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STATEMENT  OF  PROBLEM  STUDIED 


The  problem  studied  was  motivated  by  the  technological  need  to  efficiently  predict 
unsteady  transonic  flow  with  vortex  effects.  Although  much  progress  has  been  made  in  the 
last  decade1  in  the  numerical  simulation  of  transonic  flow  typical  of  that  occurring  over 
wings  of  commercial  aircraft,  these  flows  usually  have  negligible  vortical  effects.  In 
contrast,  missiles  flying  at  transonic  speeds  experience  important  effects  due  to  separation 
and,  thus,  vortical  effects.  Consequently,  the  purpose  of  the  present  research  was  to  derive 
a  reduced  set  of  governing  equations  which  properly  account  for  vortical  effects  in  unsteady 
transonic  flow,  and  to  explore  various  implementations  of  this  theory  using  well-tested 
potential  flow  solvers  as  a  basis.  The  result  will  be  efficient  and  cost  effective  computations 
for  maneuvering  missiles  and  aeroelastic  effects. 


SUMMARY  OF  MOST  IMPORTANT  RESULTS 


A  theory  was  developed  to  treat  flow  separation  and  related  vortex  effects  in  unsteady- 
transonic  flow  around  slender  bodies.2'3  This  theory  involves  the  simultaneous  solution  of  a 
potential  equation  (or  modified  Transonic  Small  Disturbance  Equation),  a  kinematic  vector 
potential  equation,  and  a  three-dimensional  unsteady  vorticity  transport  equation  for  the 
streamwise  (dominant)  component  of  vorticity.  The  theory  constitutes  a  subset  of  the  Euler 
equations,  where  flow  separation  is  modeled  by  introducing  vorticity  in  the  flow  field  by 
means  of  normal  vorticity  jets  placed  along  the  separation  line.  In  this  model,  the  location 
and  strength  of  the  separating  vorticity  are  determined  from  empirical  criteria.4  It  was 
shown  that  the  theory  can  be  regarded  as  a  unifying  theory  for  all  speed  ranges  since  it 
provides  a  framework  for  incorporating  vorticity  into  the  classic  potential  equation  for 
incompressible  as  well  as  compressible  flow.  As  a  further  refinement,  it  was  also  shown  that 
a  second  order  correction  to  the  pressure,  namely  the  inclusion  of  the  rotational  streamwise 
velocity  component,  was  critical  in  order  to  properly  account  for  the  tilting  of  lee-side 
vortices  away  from  the  body  axis. 
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An  extension  of  the  analysis5  that  included  two  components  of  vorticity  showed  that  in 
the  presence  of  constant  shed  vorticity  the  vortex  must  expand  in  size  as  it  progresses  down 
the  body.  This  extended  analysis  also  explained  the  creation  of  entropy  when  vorticity  is 
introduced  at  a  non-normal  angle  to  the  body,  thus  resulting  in  a  loss  of  total  pressure 
within  the  rolled-up  vortex.  An  important  result  is  that,  for  slender  bodies,  the  three 
components  of  vorticity  are  separable  in  magnitude. 

The  numerical  implementation  of  the  steady  version  of  the  theory  was  performed 
using  a  considerably  developed  form  of  the  TWING  potential  code.5  Results  for  bodies  at 
high  angles  of  attack  and  at  subsonic  and  transonic  Mach  numbers  were  obtained  and 
cpmpared  to  experiment  and,  for  subsonic  flow,  to  a  discrete  vortex  method.  These 
demonstrated  the  need  for  the  aforementioned  second  order  correction  on  the  pressure.  It 
was  also  found  that  numerical  errors  can  dissipate  an  isolated  vortex  by  a  considerable 
amount,  suggesting  the  need  for  further  refinement  of  the  vorticity  transport  algorithm. 

Similar  problems  were  observed  with  the  approximate  factorization  algorithm  used  in 
the  unsteady  implementation  of  the  theory.3  The  unsteady  code  was  based  on  the  well- 
tested  unsteady  transonic  potential  flow  solver  CAP-TSD.  This  code  is  capable  of  predicting 
the  unsteady  transonic  and  supersonic  flow  over  a  complete  aircraft  configuration  including 
fins,  stores,  and  pylons.7  Improvements  in  the  representation  of  the  fuselage  on  the  one 
hand,  and  the  development  of  an  appropriate  transfer  of  boundary  conditions  between  the 
true  and  computational  surfaces  on  the  other,  were  found  to  be  key  considerations  in  being 
able  to  resolve  vorticity  transport  close  to  the  body.8  Time-accurate  calculations  using  Ihis 
substantially  enhanced  version  of  the  CAP-TSD  code  were  performed  for  complete  missile 
configurations  in  the  subsonic,  transonic,  and  supersonic  flow  regimes.  The  results 
demonstrated  that  the  flow  around  realistic  angle  of  attack  configurations  could  be 
calculated  for  unsteady  transonic  flow  with  separation,  thus  showing  considerable  potential 
for  aeroelastic  computations  and  unsteady  aerodynamics. 

Due  to  the  unavailability  of  experimental  data  for  unsteady  separated  transonic  flow 
over  three-dimensional  slender  bodies,  comparisons  of  time-accurate  solutions  with 
experiment  could  only  be  performed  for  steady  angle  of  attack.  In  the  unsteady  cases,  a 
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quasi-steady  implementation  of  the  modified  Stratford9  criterion  was  used  to  pt  .diet  the 
instantaneous  separation  line.  The  question  as  to  whether  additional  time  lags  may 
originate  from  the  separation  process  itself  was  addressed  by  means  of  two-dimensional 
unsteady  Navier-Stokes  calculations  of  flow  around  a  circular  cylinder.  The  results  of  this 
effort  indicated  that  the  time  dependent  behavior  of  integral  quantities  (such  as  vorticity 
flux  and  drag  coefficient)  and  to  a  lesser  extent  of  the  separation  angle  could  be  predicted 
with  great  accuracy  using  indicia!  theory.10  The  location  of  the  separation  point  was  found 
to  be  satisfactorily  predicted  within  the  range  of  frequencies  corresponding  to  missile 
flutter.  In  particular,  it  was  shown  that  substantial  phase  delays  and  overshoots  of 
fluctuations  in  the  vorticity  flux  could  be  attained  within  the  range  of  frequencies  of  interest. 
Although  these  results  were  obtained  for  low  Reynolds  number  laminar  flow,  they  establish 
for  the  first  time  that  several  key  aspects  characterizing  the  time-dependent  behavior  of 
two-dimensional  boundary  layer  separation  may  be  predicted  using  indicial  theory. 
Whether  such  ideas  are  applicable  to  high  Reynolds  number  turbulent  flow  for  the 
development  of  fully  unsteady  separation  criteria  remains  to  be  established. 

The  unsteady  computer  code  that  was  developed  under  the  present  contract  is 
currently  being  used  for  sea-skimming  missile  applications  (NSWC  Contract  No.  N60921- 
90-C-0134),  in  which  missiles  flying  close  to  the  sea  surface  experience  periodic  gusts  due  to 
the  waves.  These  gusts  can  excite  structural  (aeroelastic)  modes  which  result  in  a 
deterioration  of  the  control  systems'  performance.11  The  unsteady  aeroelastic  effects  make 
the  code  ideally  suited  for  this  application.  The  code  is  currently  being  enhanced  to 
calculate  the  transonic  and  supersonic  flow  over  a  swept-finned  missile  of  arbitrary  roll 
angle  configuration.  This  additional  enhancement  will  be  supplemented  with  a  more 
complete  version  of  the  Transonic  Small  Disturbance  Equation  which  includes  the  necessary 
nonlinearities  for  swept  horizontal  as  well  as  vertical  shocks. 
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The  Occurrence  of  Multiple  Solutions 
for  the  "TSD-Euler"  Equation 

by 

David  Nixon 

Nielsen  Engineering  &  Research,  Inc. 
Mountain  View,  California 
94043-2287  USA 


Summary 

In  a  recent  paper  the  TSD-  Euler  equation  for  transonic  flow  has  been 
derived.  This  equation  is  similar  in  some  respects  to  the  TSD  equation  but  has 
both  entropy  and  vorticity  terms  retained.  In  this  paper  the  existence  of  multiple 
solutions  for  the  TSD-Euler  equation  is  examined  and  it  is  found  that  such 
solutions  excist  for  a  small  range  of  Mach  numbers  and  airfoil  thickness.  It  is 
found  also  that  the  addition  of  a  vorticity  flux  on  the  airfoil  surface  can  enhance 
the  appearance  of  multiple  solutions. 


"’iitroduction 

The  existence  of  multiple  or  "phantom"  solutions  in  potential  simulations  of 
transonic  flow  was  first  described  by  Steinhoff  &  Jameson  [1].  These  solutions  give 
lift  for  a  symmetric  airfoil  at  zero  degrees  angle  of  attack  and  appeared  for  a  small 
range  of  Mach  number.  These  phantom  solutions  appear  when  some  degree  of 
asymmetry,  either  in  the  algorithm  or  in  the  initial  conditions,  is  introduced  into 
the  solution  process  for  a  certain  range  of  Mach  numbers.  The  asymmetry  may 
also  be  introduced  by  a  physical  angle  of  attack.  It  appears  that  there  is  a  certain 
range  of  Mach  number  in  which  the  phantom  solution  is  preferred  to  the 
conventional  solution.  Most  of  the  early  work  on  the  occurrence  of  phantom 
solutions  concerned  steady  flow,  although  Williams  et  al  [2]  solved  the  unsteady 
TSD  equation  to  give  phantom  solutions  for  unsteady  flows.  Williams  et  al  [2l 
also  modeled  viscous  effects  by  &  boundary  layer  displacement  thickness  and  found 
phantom  solutions  for  a  slightly  different  Mach  number  range  than  for  inviscid  flow 
models.  Nixon  [3]  attempted  to  explain  mathematically  the  occurrence  of  phantom 
solutions  and  concluded  that  the  flow  must  be  transonic.  Phantom  solutions  have 
also  been  found  for  three-dimensional  flows  [4].  These  phantom  solutions  have, 
until  recently,  not  been  observed  in  equation  sets  other  than  the  potential  equation; 
although,  it  has  been  speculated  by  Williams  et  al  [2]  and  Nixon  [5]  that  they 
would  exist  for  the  Euler  equations.  All  of  these  studies  are  concerned  with 
potential  flow,  although  Salas  et  al  [6]  did  try  to  find  phantom  solutions  for  the 
Euler  equations  with  no  success.  This  paper  is  concerned  with  the  possible 
occurrence  of  phantom  solutions  in  the  Euler  equations  although  the  analysis  uses  a 
subset  of  the  Euler  equations,  called  the  TSD-Euler  equation  derived  by  Nixon  [5]. 
It  is  concluded  that  for  very  thin  airfoils  at  a  freestream  Mach  number  close  to 
unity  phantom  solutions  to  the  TSD-Euler  equation  exist.  It  is  found  also  that  the 
addition  of  vorticity  at  the  airfoil  boundary  enhances  the  probability  of  obtaining  a 
phantom  solution. 

Examples  of  Multiple  Solutions 

An  example  of  the  pressure  distribution  around  a  symmetric  Joukowski  airfoil 
at  zero  degrees  angle  of  attack  and  M»  =  0.85  is  shown  in  Figure  1.  This 
solution  was  produced  by  solving  the  transonic  small  disturbance  (TSD)  equation 
using  the  conservative  Murman-Cole  algorithm.  The  asymmetry  was  introduced  as 
an  initial  condition.  The  phenomena  appears  for  conventional  lifting  airfoils.  An 
example  [6]  of  the  variation  of  lift  coefficient,  Cl,  with  angle  of  attack  a,  for  a 
RAE2822  airfoil  is  shown  in  Figure  2.  It  may  be  observed  that  at  a  =  0.6 
degrees,  Cl  can  have  a  value  of  about  0.5  or  about  1.4.  It  is  interesting  to  note 
that  if  the  phantom  solution  is  realizable  the  value  of  Cl  would  be  three  times  the 
"conventional"  value. 

TSD-Euler  Equation 


One  possible  reason  for  the  inability  to  find  phantom  solutions  for  the  Euler 
equations  is  that  the  algorithms  for  solving  the  Euler  equations  do  not  satisfy  these 
equations  on  the  airfoil  surface  since  the  tangential  velocity  on  the  surface  are  found 
by  extrapolation  from  the  interior.  It  is  noted  by  Nixon  [3]  that  it  is  the  behavior 
of  the  potential  equation  on  the  surface  that  gives  rise  to  the  eigensolution 
necessary  to  obtain  phantom  solutions.  To  examine  this  aspect  further,  Nixon  [5] 
derived  a  "small  disturbance"  version  of  the  Euler  equations,  the  "TSD-Euler" 


equation,  which  includes  the  first  order  effects  of  entropy  and  vorticity.  This 
equation  reduces  to  the  classic  transonic  small  disturbance  (TSD)  equation  (for 
which  phantom  solutions  can  be  computed)  in  the  absence  of  vorticity  and  has  the 
advantage  that  the  basic  equation  is  satisfied  everywhere,  thus  allowing  the 
emergence  of  the  necessary  eigensolution. 

The  TSD-Euler  equation  [5]  is  given  by 

[A-ku?]  +  v  =  (  §  ) 

1  2  Jx  7  (1) 

where  u  and  v  are  perturbation  velocities  in  the  x  and  y  directions,  respectively, 
nondimensionalized  with  respect  to  freestream  values.  S  is  the  entropy  change,  R  is 
the  gas  constant  and 


P2  =  1  -  ul‘,  k  -  (7  +  1)  M 2 


(2) 


where  Moo  is  the  freestream  Mach  number, 
and  f  is  a  vector  potential,  then 

u  =  + 

V  =  - 


If  ^  is  a  perturbation  velocity  potential, 


\ 


and 


(4) 


where  the  vorticity  w  is  related  to  the  entropy  S  by  Crocco's  theorem.  The 
entropy  generated  at  a  shock  is  given  by 


(«*2-i)3  (5) 

R  (7  +  1) 

where  M+  is  the  Mach  number  just  ahead  of  the  shock  and  is  related  to  u.  The 
entropy  generated  by  the  shock  is  constant  with  respect  to  x  aft  of  the  shock.  The 
tangency  boundary  condition  is  represented  by  the  thin  airfoil  approximation  so  that 
the  equation  and  its  boundary  condition*;  reduce  to  the  classic  TSD  formulation  in 
the  absence  of  vorticity.  The  tangency  boundary  condition  is  given  by 

-  iy  =  y>) 

y  =  40  t  (6) 

where  y  =  ys(x)  denotes  the  airfoil. 
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Equation  (1)  with  its  boundary  condition,  Equation  (6),  and  the  vorticity 
equation,  can  be  written  in  similarity  form  to  give 


(«  -  id  )  +  =  ■  27  a  k  P2  aS  K1  "  “  +)3] 


2  Jx  y  (7  +  1)‘ 
v(x,*o)  =  K  ys'(x) 


(7) 

(8) 


u  =  K  +  P2  ?-  5  *  =  ?-  "  t 


(9) 


/>2  y__  ♦  yxx  =  - » =  -21-  ^  ~= [(1_u  +)3] 

yy  (7+i)  3y 


(10) 


where 


^  f  W  =  k/^3  uf  y  =  fiy  and  u+  is  the  value  of 


u  just  ahead  of  the  shock.  In  the  derivation  of  Equation  (10)  the  first  order 
approximation  to  Crocco‘s  Equation  has  been  used;  thus 


W 


i  _a  (s/r) 
r7m»2  9y 


(U) 


If  ry  =  y  where  T  is  the  thickness  of  the  airfoil  then  a  transonic  similarity 
s  s 

parameter,  K,  can  be  defined  as 


K  =  kr/^3  (12) 

The  transonic  similarity  parameter,  K,  is  constant  for  combinations  of  Mach 
number  and  thickness;  if  the  airfoil  is  very  thin  the  Mach  number  approaches  unity 
and  the  entropy  production  vanishes  as  r+o  and  /?2+o.  Thus  for  a  given  K  it  is 
possible  to  approach  the  TSD  equation  arbitrarily  closely  by  varying  the  thickness 
and  the  Mach  number. 

The  set  of  equations,  Equations  (7-10),  can  be  solved  using  the  same  algorithm 
and  treatment  of  the  boundary  conditions  used  in  the  TSD  example,  thus  ensuring 
that  a  phantom  solution  can  be  obtained  as  /7+0. 


Addition  of  Vorticity 

If  any  vorticity  is  introduced  into  the  two  dimensional  flow  the  vorticity  vector 
will  be  aligned  with  the  spanwise  direction  and  the  magnitude  of  the  vorticity  will 
not  change  in  this  direction.  According  to  the  equations  of  Klopfer  and  Nixon  (7) 
this  vorticity  is  convected  like 


(mu)x  +  (wv)y  =  0  (13) 

to  a  fi:st  approximation.  The  velocity  induced  by  the  vorticity  is  given  by 
Equation  (4)  or  Equation  (10).  In  these  models  vorticity  is  defined  as 

vx  “  uy  =  "  (14) 

Using  Equation  (14)  and  Equation  (13)  gives  an  alternative  representation  of  the 
vorticity  transport,  namely, 


»« *  v  -  *2  *  •/  -  v  (15) 

where 

Ax  =  w  YJ  Ay  =  -  «  u  (16) 

A  suitable  boundary  condition  for  Equation  (15)  is  that  normal  derivatives  of  A  are 
specified  on  the  boundaries,  that  is  Neumann  boundary  conditions  are  used. 

Equation  (15)  can  be  solved  using  a  simple  central  differencing  scheme  with 
the  vorticity  terms  on  the  right  hand  side  set  retarded  by  one  iteration.  Equation 
(15)  can  be  solved  iteratively  with  Equations  (l),  (3),  (4),  and  (5). 

In  one  of  the  examples  discussed  in  the  next  section  a  small  amount  of 
vorticity,  comparable  in  magnitude  to  the  shock  generated  vorticity,  is  introduced 
over  a  small  part  of  the  chord  (usually  over  three  grid  points).  The  amount  of 
vorticity,  w,  is  symmetric  on  both  surfaces  of  the  airfoil;  the  symmetry  or 
asymmetry  of  u  depends  on  whether  a  phantom  solution  is  present. 


Results  from  the  TSD-Euler  Equation 

It  has  been  found  that  phantom  solutions  for  the  TSD-Euler  equation, 
Equation  (l),  for  the  NACA00XX  series  airfoil,  can  only  be  obtained  for  a  very 
thin  airfoil  at  a  freestream  Mach  number  close  to  unity.  This  is  in  keeping  with  a 
suggestion  of  Williams  et  al  [2].  An  example  for  a  0.00076%  thick  airfoil  at 
M»  =  0.9945  is  shown  in  Figure  3.  The  lowest  Mach  number  for  which  a 
phantom  solution  could  be  obtained  is  0.986  for  an  airfoil  of  thickness  0.003%. 
These  results  indicate  that  the  occurrence  of  the  phantom  solutions  in  the  Euler 
equations  is  not  likely  to  be  of  practical  interest. 
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As  an  experiment  a  vorticity  source,  equivalent  to  a  heat  sink,  was  placed 
upstream  of  the  >  shock  waves  on  the  airfoil  surface.  The  entropy  induced  by  this 
source  is  approximately  the  same  magnitude  of  that  produced  by  the  shock  but 
with  the  opposite  sign;  (S/R)  is  of  the  order  of  10*4.  Using  this  device  it  proved 
possible  to  generate  a  lifting  solution  for  a  Mach  number  of  0.945  for  a  2  1/2% 
thick  airfoil;  the  solution  is  shown  in  Figure  4.  This  is  a  much  more  interesting 
proposition  than  the  previous  cases.  It  can  be  inferred  from  this  result  that  a 
small  amount  of  vorticity  can  greatly  increase  the  chance  of  getting  the  phantom 
solution. 

If  the  phantom  solutions  found  for  the  TSD-Euler  equations  could  be  made 
physically  realizable  they  are  only  practical  if  a  vorticity  source  is  inserted  into  the 
airfoil  boundary,  indicating  that  a  control  device  is  necessary  to  obtain  these 
solutions.  The  fact  that  phantom  solutions  have  been  found  for  flows  with  shock- 
induced  and  distributed  vorticity  suggests  that  the  solutions  can  be  obtained  for  the 
Navier-Stokes  equations.  The  calculations  reported  by  Williams  et  al  [2]  for 
boundary  layer  flows  reinforce  this  point.  It  is  suggested  that  the  phantom  solution 
are  realizable  and  that  these  could  provide  radically  increased  aerodynamic  efficiency 
of  conventional  flows  as  may  be  observed  from  the  result  shown  in  Figure  2. 

Concluding  Remarks 

An  equation  with  many  similarities  to  the  Euler  equations  has  been  examined 
for  the  occurrency  of  multiple  or  phantom  solutions.  It  is  concluded  that  although 
these  solutions  do  exist  it  is  for  a  much  more  restrictive  range  of  Mach  numbers 
than  the  potential  equation.  If  a  small  amount  of  vorticity  is  introduced  into  the 
flow  through  the  airfoil  surface,  upstream  of  the  shock,  then  phantom  solutions  can 
be  found  for  a  much  larger  range  of  Mach  number.  Since  this  model  has  vorticity, 
entropy  and  shock  waves  represented  it  is  speculated  that  phantom  solutions  could 
be  found  for  the  Navier-Stokes  equations.  The  fact  that  such  solutions  of  the 
potential  equation  coupled  with  a  boundary  layer  model  serves  to  reinforce  this 
point. 


If  the  Navier-Stokes  equations  have  phantom  solutions  similar  to  those  found 
for  the  potential  or  TSD-Euler  equations  then  it  may  be  possible  to  exploit  this 
"new"  aerodynamic  phenomena  to  radically  improve  aerodynamic  efficiency.  In  fact 
there  may  be  many  instances  of  phantom  solutions  involving  more  than  one  location 
of  concentrated  vortices,  shock  waves,  or  shear  surfaces  that  could  be  used  to 
revolutionize  aircraft  flight  if  these  phenomena  can  be  made  to  occur  and  stabilized 
in  a  real  flight. 
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Introduction 

In  recent  years  it  has  been  possible  to  predict  the 
unsteady  transonic  flow  around  a  wing,  especially  those 
typical  of  commercial  aircraft,  in  a  fairly  efficient 
manner.  Frequently,  the  computer  codes  that  are  used 
are  based  on  potential  theory  and  are  considerably  faster 
than  a  corresponding  calculation  using  the  Euler 
equations.  However,  these  methods  have  been 
developed  for  aircraft  and  are  not  really  applicable  to 
missile  flows  where  vorticity  effects  are  important;  the 
potential  approximation  cannot  predict  the  effects  of 
vorticity  other  than  by  representing  a  vortex  wake  by  an 
infinitely  thin  sheet  which  is  excluded  from  the 
computational  domain.  This  model  is  complicated  and 
may  not  be  a  viable  option  for  routine  calculations 
around  real  aircraft  or  missile  configurations  because  the 
geometry  of  the  vortex  sheet  can  get  quite  complex. 

Missiles  generally  have  low  aspect  ratio  fins  and 
slender  bodies,  neither  of  which  are  features  of 
commercial  aircraft;  this  may  cause  problems  in  grid 
generation.  The  airflow  over  missiles  is  dominated  by 
vortex  effects;  this  is  in  contrast  to  the  attached  flow  over 
a  commercial  aircraft.  For  steady  missile  flow  a  variety 
of  prediction  methods  are  available,  ranging  from  panel 
methods,  with  the  addition  of  nonlinear  vortex 
dynamics,  to  the  Euler  equations  or  Navier-Stokes 
equations.  The  Navier-Stokes  equations  will  model  flow 
separation  and  track  the  resulting  vorticity  but  are  only 
as  accurate  as  the  turbulence  model  used  in  the 
calculation.  This  is  critical  in  defining  the  separation 
line.  The  prediction  methods  based  on  the  Euler 
equations  need  some  empirical  criteria  to  predict -the 
location  and  strength  of  the  separating  vortex  sheet;  the 
Euler  equations  will  then  model  the  transport  of  this 
vorticity  reasonably  accurately.  A  detailed  discussion  of 
the  methods  used  to  predict  steady  flows  over  missiles  is 
given  in  Reference  (1).  In  order  to  predict  unsteady  flow 
over  a  missile  to  adequate  accuracy,  the  method  must  be 
able  therefore,  to  treat  unsteady  separated  flow. 
Furthermore,  if  it  is  necessary  to  predict  the  steady  flow 
at  transonic,  rather  than  supersonic,  speeds  the  method 
mi  st  be  computationally  efficient  since,  unlike 
^upeisuiiiv.  fiuVv  Which  can  usua  marching  piuccdutu, 
transonic  flow  methods  require  numerous  sweeps  of  the 
computational  domain  leading  to  large  computational 
times. 
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It  would  be  ideal  to  use  the  Navier-Stokes 
equations  to  model  unsteady  transonic  flow  around 
missiles,  but  there  are  several  difficulties  with  this 
approach.  The  most  obvious  difficulty  is  the  computer 
time  required,  which  is  several  orders  of  magnitude 
greater  than  that  required  for  a  potential  calculation. 
Even  if  a  dedicated  supercomputer  were  available  for 
such  calculations,  the  computation  of  the  necessary  flow 
separation  might  be  inaccurate  because  of  the  inherent 
inaccuracy  of  many  present  turbulence  models.  The  next 
best  approach  is  to  use  the  time  dependent  Euler 
equations  which  require  some  empirical  factor  to  initiate 
separation  and  compute  the  shed  vorticity.  This 
empiricism  would  be  a  development  of  the  steady  flow 
criteria  described  in  Reference  (1).  If  the  separation  line 
and  shed  vorticity  can  be  predicted,  then  the  Euler 
equations  would  be  a  viable  model.  However,  since 
missiles  are  slender  some  further  approximations  can  be 
made.  This  paper  is  concerned  with  developing  a 
method  of  predicting  the  pressure  distribution  on 
missiles  at  transonic  and  supersonic  speeds.  The 
concepts  developed  in  this  work  will  be  transferred  to  an 
unsteady  analysis.  The  approach  is  to  make  as  much  use 
of  the  existing  technology  as  possible.  The  final  goal  is  to 
develop  a  computer  code  capable  of  predicting  the 
unsteady  transonic  flow  about  missiles  for  use  in 
aeroelastic  calculations  or  during  maneuver.  The 
present  paper  is  an  account  of  the  first  phase  of  the  work, 
namely  the  development  of  a  steady  flow  variant  of  the 
method. 

The  approach  is  based  on  experience  gained  in 
steady  subsonic  and  supersonic  flow  predictions  for 
.missiles,  in  particular,  the  fact  that  if  a  separation  line 
and  the  strength  of  the  vorticity  introduced  at  that  line 
can  be  estimated  empirically,  the  governing  equations 
(such  as  the  Euler  equations)  will  repiesent  the  transport 
of  this  vorticity  accurately.  Because  the  computational 
time  for  the  Euler  equations  for  an  unsteady  calculation 
is  considerably  greater  than  that  for  a  potential 
calculation,  a  simplified  model  is  used. 

The  basic  equations  ere  those  derived  by  Klopfer 
and  Nixon2  and  are  a  subset  of  the  Euler  equations.  For 
slender  bodies,  the  five  Euler  equations  are  reduced  to 
an  equation  similar  to  the  three-dimensional  unsteady 
potential  equation,  a  vorticity  transport  equation  and  a 
two-dimensional  Poisson  equation.  Apart  from  the 
reduced  set  of  equations,  a  second  advantage  is  that  one 
of  the  equations  is  almost  identical  to  the  potential 
equation  for  which  there  are  well  tested  computer  codes. 
This  allows  the  development  of  a  prediction  method 
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based  considerably  on  proved  technology.  A  derivation 
of  the  equations  is  given  in  Reference  2. 

Basic  Equations 

From  Reference  2,  the  basic  equations  are  as 
follows: 

Vorticity  Transport 

u»nx  +  <vnV  +  (wn)z  85  0  0) 

Rotational  Crossflow 

Ayy  +  ^“zz  =  (2) 

Mass  Conservation 

(pu)x  +(pu)y  +(pw)z  =0  (3) 

where 

u  =  tfx;v  =  *y  +  Az; w  =  $z  -Ay  (4) 

and 


f  (7“1)  M 
.1**-  * 


2  2  2. 
(1“U  -v  -w  ) 


1 

7-1 


(5) 


4>  is  the  velocity  potential,  A  is  the  vector  potential,  p  is 
the  density,  Mw  is  the  freestream  Mach  number,  0  is  the 
vorticity  in  the  streamwise  direction,  and  u,  v,  and  w  are 
velocities  in  the  coordinate  directions  x,  y,  z  respectively. 

Governing  Equations  in  Transformed 
Curvilinear  Coordinates 


derived.  The  total  velocity,  V,  is  defined  as  the  sum  of 
potential  and  rotational  velocity  components,  p  and 
VR,  respectively: 


V^’+VR  (7) 

where 

=v*  (8) 

and 

VR  =  VxA  (9) 

with 

X  =  (A,0,0)t  (10) 

In  Cartesian  coordinates: 

Vp  =(tx.4y,  tz)r  (Ha) 

and 

\>R  =  (0,  Az,  -  Ay)T  (lib) 

Note  that 

Vxtfp  =0  (12a) 

and 

7.tfR=0  (12b) 

In  other  words,  the  potential  velocity  field  is  irrotational 
and  the  rotational  velocity  field  is  divergence  free. 

The  continuity  equation,  (3),  in  vector  form  is 


V  .  (pV)=0 


(13) 


In  order  to  treat  complex  geometries,  the 
governing  equations  are  transformed  to  a  curvilinear 
body  conforming  coordinate  system  which  is  orthogonal 
in  the  crossflow  plane: 


Combining  Equations  (7)  -  (13)  gives 
V.  (p7tf)  =  S 

where 


S  =  -tfR  ♦  Vp 


(14) 

(15) 


y=y(*/T},$)  (6) 

z  =  z(£,tj,$) 

Note  that  in  the  transformed  coordinates,  Pt  is  the 
streamwise  direction,  the  ij  and  $  directions  are  a  bai/is 
for  the  crossflow  plane,  and  Vrj  •  7$  =  0.  The  rj  and  $ 
directions  are  parallel  and  normal,  respectively,  to  the 
body  surface. 

Before  presenting  the  transformed  equations,  an 
alternate  formulation  of  the  continuity  equation  will  be 


Equation  (14)  is  the  full  potential  equation  with  a  source 
term.  It  can  be  transformed  from  Cartesian  coordinates 
to  strong  conservation  law  (SCL)  form.  An  existing  code 
is  used  to  solve  Equation  (14)  in  SCL  form,  and  will  be 
described  below. 

The  vorticity  transport  equation  transforms  to  the 
following  SCL  form: 

(16) 

where 


A  ft 

n  -  2  (i7) 

and  the  contravariant  velocities  are  defined 

V  =  r}x  +  ijy  v/ua>+  n2  w/uM  (18a) 

w  =  $x  +  $y  w/ua+  Jz  w/uw  (18b) 

Formulas  for  the  Jacobian,  J,  and1  the  metric  terms  rjx , 
ijy ,  $x,  etc.  can  be  found  in  Reference  3  for  a 
generalized  curvilinear  transformation.  The 
transformation,  Equation  (6),  permits  simplification  of 
some  of  these  terms. 

The  vector  potential  equation,  (2),  transforms  to: 

aAT?T?  +  BA$$  +  7Ai}+  6A$  =  ‘n  <19) 

where 

«~r}y  +v\ 

7  =  7>yy  +  7*z z 


Since  the  transformation  is  orthogonal  in  the  rj-$  plane, 
the  cross-derivative  term,  does  not  appear. 

To  summarize,  the  governing  equations, 
continuity,  vorticity  transport,  and  the  vector  potential 
equations  have  been  transformed  to  a  curvilinear 
coordinate  system  and  are  given  by  Equations  (14),  (16), 
and  (19),  respectively.  The  following  section  discusses 
boundary  conditions  used  for  these  equations. 


Boundary  Conditions 

Boundary  conditions  for  the  velocity  potential,  <j>, 
the  vorticity,  ft,  and  the  vector  potential.  A,  are  required 
on  the  body  surface  and  at  the  far  field  boundary.  All 
boundary  conditions  will  be  described  next. 

Body  Surface 

Away  from  the  separation  locations,  flow  tangency 
is  applied  on  the  body  surface: 


where  )  is  a  normal  velocity  component.  Equation 
(21a)  indicates  that  either  or  V*  can  be  arbitrarily 
specified.  We  choose 

V£=0  (21b) 

since  this  condition  is  built  into  the  existing  full  potential 
code.  With  VjJ  chosen,  Equation  (2la)  requires  that 

V*=0  (22) 

In  a  locally  orthonormal  coordinate  system  in  the 
crossflow  plane  with  n  and  t  (the  unit  tangential  to  the 
surface)  as  the  basis  vectors,  the  normal  component  of 
the  rotational  velocity  becomes 

VR  =-|i-  (23) 

and  the  tangential  component  is: 

VR  =  &  (24) 

Equations  (22)  and  (23)  provide  the  body  boundary 
condition  for  the  vector  potential: 

$-=0  (25) 

or 

A  =  A0  =  const.  (26) 

along  the  body  (away  from  the  separation  point). 

Vorticity  is  introduced  into  the  flowfield  only  in 
the  vicinity  of  the  separation  locations  (see  discussion 
below).  Consequently,  away  from  the  separation  points 

ftl  -  0  (27) 

body 

is  applied. 

Far-Field 

In  the  farfield,  freestream  conditions  prevail, 

V  =  Va  (28) 

or  in  component  form: 

4>  I  ~  u  (29a) 

X  I  co 


V  •  n  I  -  0 
'body 


(20) 


rli  - 

3n  dfc  » 


v 

n,  oo 


(29b) 


where  It  is  the  outward  unit  normal  to  the  surface.  This 
requires  that 


Vj  +VR  =0  (2la) 


d<(>  3A* 

[■*£  +  -r-1  -  V 

o t  on  co  t,» 


(29c) 
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Where  $^-and  |^-are  local  normal  and  tangential 
derivatives  with  respect  to  the  freestream  boundary 
surface. 


Note  that  either  a  Neuman  condition  on  A,  which 
specifies  or  a  Dirichlet  condition,  which  indirectly 
specifies  can  be  applied  at  the  far  field  boundary, 
but  not  both.  Having  chosen  either  of  these  conditions, 
the  velocity  potential  boundary  condition  must  then  be 
modified  to  ensure  satisfaction  of  Equation  (29).  A  third 
alternative  is  chosen.  The  far  field  boundary  is  placed 
far  enough  away  from  the  body  such  that 


VR=|£-«0  (30) 


Then  by  specifying  a  Dirichlet  condition 


f  V  a  dl  -  ¥§  X  (36) 

J  n  2 

a 

is  satisfied,  where  the  integration  is  taken  along  the 
surtace  coordinate,  1,  in  the  vicinity  of  the  separation 
location,  a  <  1  <  b.  Note  that  an  empirical  factor,  X,  has 
been  included  to  reduce  the  boundary  layer  vorticity 
flux.  This  "vortex  reduction  factor"  is  widely  used  in 
discrete  vortex  methods  to  provide  better  agreement 
between  predicted  and  measured  values.  The 
recommended  value4  lies  in  the  range  0.6  £  X  <  1.0. 


The  separation  line  is  determined4  by  the  Stratford 
criterion  which  for  turbulent  boundary  layers  states  that 
separation  occurs  when 


C 

P 


x  10 


-0.1 


0.35  aina  (37) 


which  implies 

V*=|*.  =  0  (32) 

The  boundary  condition  for  the  velocity  potential  thus 
becomes 

V  (33) 


where  Cp  is  the  pressure  coefficient,  £  is  the  length  of 
the  boundary  layer  in  the  streamwise  direction,  and  Re  ^ 
is  the  Reynolds  number  based  on  this  length.  The 
separation  criteria  and  the  strength  of  the  shed  vorticity 
have  been  used  successfully  in  many  subsonic 
examples  4  As  a  first  step  the  same  criteria  are  used 
here. 


and  this  is  already  implemented  in  the  existing  full 
potential  code.  Equations  (31)  -  (33)  indicate  that  the  far 
field  flow  is  irrotational,  consequently: 

O  |  -  0  (34) 

00 

is  also  specified  at  the  far  field  boundary. 

Separation  condition 

The  vorticity  flux  passing  through  a  plane  normal 
to  a  two-dimensional  boundary  layer  is  given  to  the 
order  of  approximation  of  boundary  layer  theory  by 

|vj)dy-|v|^dy-^  (35) 

where  Ve  is  the  velocity  at  the  edge  of  the  boundary 
layer,  5. 

At  the  separation  point,  it  is  assumed  that  this 
vorticity  is  "injected"  into  the  inviscid  flow  field. 
Consequently,  a  vorticity  flux  boundary  condition  is 
applied  to  the  vorticity  transport  equation  at  the 
separation  location.  The  flux  (Vn  ft)  is  applied  at  the 
body  surface  such  that 


It  is  shown  in  Reference  2  that  the  present 
formulation  is  equivalent  at  subsonic  speeds  to  the 
"vortex  cloud"  method  used  by  Mendenhall  and 
Perkins.4  The  main  difference  is  that  the  vorticity 
transport  equation.  Equation  (1),  is  solved  by 
representing  the  vorticity  flux  by  discrete  vortices  which 
are  transported  with  the  flow.  The  velocity  induced  by 
the  vortices  is  computed  from  the  Biot-Savart  law.  The 
strength  of  the  vortices  introduced  at  the  separation  line 
is  determined  by  the  same  criteria  as  the  present  work, 
that  is  Equation  (36).  However,  the  location  of  a  vortex 
that  has  just  been  introduced  must  be  determined,  which 
is  achieved  by  ensuring  that  the  separation  line  is  a 
crossflow  stagnation  point.  Since  the  discrete  vortex  is 
the  integral  of  the  vorticity  over  an  element,  the  location 
may  be  regarded  as  the  centroid  of  the  vorticity  in  the 
element.  If  the  vorticity  flux  in  a  time  interval  is 
uniform,  then  the  centroid  is  at  the  center  of  the  element; 
there  is  no  flexibility  in  changing  the  centroid  other  than 
by  changing  the  rate  of  change  of  vorticity  flux  in  each 
time  interval  which  seems  to  be  a  nonphysical  artifice. 
Another  difference  between  the  present  method  and 
vortex  cloud  is  that  discrete  vortex  methods  do  not 
automatically  converge  as  the  number  of  vortices  is 
increased;  this  introduces  an  element  of  art  into  the 
procedure  which  cannot  be  duplicated  in  the  present 
method. 
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A  Second  Order  Correction 


In  the  work  of  Mendenhall  and  Perkins,4  it  is 
necessary  to  introduce  the  rotational  increment  in  u  due 
to  the  rotation  of  the  vortex  away  from  the  axial 
direction.  This  additional  term  is  necessary  to  obtain  any 
kind  of  adequate  results  and  is  presented  as  an  empirical 
correction.  It  is  necessary  to  include  this  correction  in 
the  present  theory  but  it  is  instructive  to  examine  the 
idea  further. 

From  the  definition  of  the  various  vorticity 
components 


wRy  -v\  =0, 

U\  -wRx  =«2 

VRX  -uR  =03 


(38a) 

(38b) 

(38c) 


where  the  superscript  "R"  denotes  a  rotational 
component,  ftj  is  the  vorticity  component  in  the  axial 
direction,  is  the  component  in  the  y  direction,  and  ft  3 
is  the  component  in  the  2  direction.  In  the  computations, 
fij  is  computed.  From  Equation  (38b) 

T>  f  Z  ft 

U  -  ]  (w*  +  ft2)dz  +  f(x,y)  (39) 

where  f(x,y)  is  a  function  due  to  i^e  integration  and  on 
differentiation 


R 

UY  " 


(vr  +  n  )  dz  +  f 
xy  2y  y 


(40) 


From  Equation  (38a)  and  (38c) 

R  0  ,  R 

w  -  0,  +  v 
y  1  z 


(41) 


2 

uR  -  vR  -  ft-  I  (vR  -  ft.  )dz 
y  x  3  J  xz  3z 


(42) 


R 


u 


+ 


ft,  Idz 


(44) 


From  the  analysis  of  Klopfer  and  Nixon2  it  may  be 
shown  that 


WR  -O.Ceu^) 


ft2  -a(eOj)-  0(e3) 

Herce  a  First  approximation  is 

R  (2  R. 
u  -  ]  wxdz 


(45) 


(46) 


The  point  of  this  analysis  is  that  the  "empirical 
correction"  of  Mendenhall  and  Perkins4  is  actually  a 
consistent  second  order  extension  of  the  theory.  This 
extra  component  of  velocity  is  incorporated  into  the 
pressure  relation 

v^{[l+^M»(l-q2)f/7 1-1 }  (47) 

'  CO 


where 

q2 -(up +uR)2  +  (vp+vR)2  +(wp+wR)2  (48) 

Numerical  Methods 

The  governing  equations,  continuity,  vorticity 
transport,  vector  potential  Equations  (14)  -  (19),  and  the 
Bernoulli  equation  are  a  coupled,  non-linear  system  and 
are  solved  in  an  iterative  manner.  The  numerical 
method  used  for  each  equation  and  the  overall  solution 
algorithm  for  the  system  are  described  in  this  section. 

Continuity  Equation 


and  on  substitution  of  Equations  (41)  and  (42)  into 
Equation  (40) 

f  “  (  [yR  “  8-  -  ft-  -  vR  -  ft.  I dr7  (43) 

y  j  I  xz  3z  lx  xz  2yJ  '  v 


Equation  (14)  with  S  =  0  is  the  transonic  full 
potential  equation.  An  existing  code,  TWING3 ,  which 
solves  the  full  potential  equation  in  SCL  form  for  a  v/ing 
was  modified  to  treat  a  general  body  and  to  include  the 
source  term,  S,  and  the  terms  Ay  and  Az  in  the  Bernoulli 
equation. 


TWING  uses  an  implicit,  approximate  factorization 
(AF)  shock-capturing  scheme.  Central  differencing  is 
Since  used  at  sub.sonic  points  and  an  equivalent  upwind 

density  biasing  is  applied  at  supersonic  points.  For  a  full 
fyx  +*hy  +^3z  description  of  the  numerical  method  used  in  TWING,  see 

Reference3. 

by  definition  it  follows  that 

The  AF  scheme  for  the  full  potential  equation  can 
be  written 
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NCrt  +o>Un  =  0 


(49) 


where  =  0  represents  the  full  potential  equation  in 
operator  form.  The  is  a  relaxation  parameter;  Cn  is  the 
correction  term  (^n  + 1  -  ),  L^n  is  the  residual,  and  N 

is  an  operator  representing  the  AF  splitting.  The  source 
term  was  incorporated  into  the  AF  scheme  by  replacing 
the  residual  L$n  with(L$n  -Sn)  giving 

NCn  -c o(Un  -Sn)  =  0  (50) 

and  Sn  is  evaluated  according  to  Equation  (15).  The 
velocity  terms  appearing  in  the  Bernoulli  Equation  (5) 
were  modified  to  include  Ay  and  Az  as  indicated  in 
Equation  (4).  No  further  modifications  were  required. 

Vorticity  Transport  Equation 

The  vorticity  transport  equation  is  hyperbolic  with 
{  as  the  time-like  variable.  Consequently,  a  marching- 
type  procedure  is  used  to  solve  this  equation.  Central 
differencing  is  used  with  a  Beam  and  Warming-type  5 
solution  algorithm  that  is  derived  next. 


V  >1*  " 1/2  {  (  >1,1, 


Jc+1 


'  (  *i, j,k-l} 

Equation  (52)  can  be  rearranged: 


I  +  ~  [«  ViM  +  $,Wj 


Vi+1  ri+U 


*1+1 


-  H  [\vi +  5jwi]j  ai 


Equation  (54)  is  then  factorized  giving, 


(53b) 


(54) 


(I  +  t  5/i+i> (I  +  *  6j*W  «i+i 


The  vorticity  transport  equation  (16)  is  first 
rewritten  as 


n{  -  - [  (To)  _  +  (wn>{|  (si) 

Trapezoidal  differencing  is  used  in  the  {-direction,  and 
second-order  central  differencing  in  the  rj  and  { 
directions,  resulting  in: 


fti+l  “  ni 
A{ 


-  -  1/2 


6  (Vft)  +  6,(Wfi) 
V  s 


i+1 

(52) 


+  [fi^Vft)  +  6 ^  (W^)  ] 


(55) 


(X  -  M  SVL)  (I  -  *  6^)  ai  +  Dt 


Note  that  this  equation  is  linear  in  ft,  since  provisional 
values  of  + 1  and  Wj  + 1  are  known  when  Equation 
(55)  is  solved  (see  discussion  below  concerning  the 
global  iteration  strategy).  The  right  hand  side  of 
Equation  (55)  has  also  been  factored  to  minimize  the 
factorization  error. 

Dj,  in  Equation  (55),  represents  an  added  explicit 
fourth-order  dissipation  term  which  is  necessary  to 
stabilize  the  integration  in  {.  The  dissipation  also 
provides  a  spatial  smoothing  of  the  vorticity  which 
inhibits  the  odd-even  decoupling  phenomena, 
characteristic  of  central-differenced,  convective-transport 
equations. 


where  "i"  is  the  grid  index  in  the  {-direction;  "j"  and  ”k," 
the  indices  in  the  tj  and  $  directions,  respectively  are 
implied  in  Equation  (52).  The  central  difference 
operators  are 

V  ’i*  “ 1/2 1  < 

(53a) 

(  ’i.j-l.k} 

and 


The  formula  used  for  Dj  is  taken  from  Reference  6: 
Di  =  <J;i(84,  +  s4jKnj)i  (56) 

where  54  and  6  j  are  centered,  fourth  differences,  e.g.: 

54n  =  0.2  +6^-40;.,  +  ftjt2  (57) 

'< 

and  e,  the  fourth  order  smoothing  coefficient,  is  a 
constant. 

The  fourth  order  difference  terms.  Equation  (57), 
are  modified  near  the  boundaries  and  are  discussed 
below  in  the  results  section. 
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To  advance  Equation  (55),  the  following  two  step 
procedure  is  executed. 

1:{I  +  ^5,,V|+,)  flA;+1  =Rj  (58a) 

2:  (I  +  fl65WM)  SiAit,  =S;+1  (58b) 

where  R.  =  right-hand-side  of  Equation  (55).  Steps  1  and 
2  require  the  inversion  of  scalar  tridiagonal  matrices. 

Vector  potential  equation 

The  vector  potential  is  a  Poisson  equation  that  is 
solved  separately  in  each  crossflow  plane.  Second  order 
central  differencing  is  used  for  both  first  and  second 
order  derivatives  in  Equation  (19).  Successive  over- 
relaxation  is  used  as  the  iteration  schema 

The  second  order  metric  derivatives  appearing  in 
Equation  (19)  are  evaluated  by  using  the  chain  rule,  e.g.*. 

Vyy  +(»?y)$$y  (59) 

(recall  that  £  *  f(y)).  Central  differencing  is  used  to 
approximate  the  terms  like  (rjy)^ 

Treatment  at  grid  singularities 

Two  different  mesh  topologies  were  used  in  the 
study,  an  H-grid  and  a  polar  grid.  At  the  grid 
singularities,  the  dependent  variables  were  taken  as 
simple  averages  of  neighboring  values  in  the  same 
crossflow  plane. 

Global  iteration  strategy 

The  solution  of  the  coupled  system  of  Equations 
(14),  (16),  and  (19)  is  achieved  by  sequentially  iterating 
each  individual  equation  in  turn.  A  single  pass  through 
all  the  equations  is  a  global  or  outer  iteration;  inner 
iterations  are  performed  on  the  velocity  potential  and 
vector  potential  equations.  Because  the  inner  iteration 
algorithms  were  designed  to  be  stable  and  convergent, 
the  global  iteration  process  is  expected  to  be  similarly 
well  behaved. 

One  global  iteration  consists  of  the  following  steps: 

1.  Iterate  on  the  velocity  potential  equation  (50). 
Initially  S  =  Az  =  Ay  =  0. 

2.  Determine  the  incipient  separation  plane  £  -  £{  s 
using  the  Stratford  criterion  in  Equation  (37).  This  is 
the  first  upstream  plane  where  vorticity  will  be 
shed. 

3.  March  the  vorticity  transport  and  vector  potential 
equations  downstream  from  £is  to  £max  where 
*max  is  end  value  of  £. 


*!  “  *is 

(b)  Determine  the  separation  location  from  the 
Stratford  criterion  and  evaluate  the  vorticity 
flux  to  be  "shed"  to  £.  crossflow  plane 
according  to  Equation  (36). 

(c)  Advance  vorticity  from  plane  £.  to  £i  +  1 . 
Since  vorticity  transport  and  vector  potential 
equations  are  coupled  through  VR ,  this  is  an 
iterative  procedure,  n  =  1. 

(i)  Step  vorticity  transport  Equation  (51) 
from  plane  £j  to  £l+r 

(ii)  Iterate  vector  potential  Equation  (19)  for 
£j+i  plane. 

(iii)  Calculate  total  velocity  for  £•  .  ■,  plane: 

(iv)  Check  for  convergence  on  total  velocity  in 
£id  P'ane:  if 

|  V f  +  ^  -  V1}'}  I  -  0  iteration 
converged,  otherwise  continue  iterating 
(n  =  n  + 1,  go  to  (i)  above). 

(d)  If  at  last  crossflow  plane,  £i  +  1  =  £m  a  x 
continue,  otherwise  continue  marching  vorticity 
downstream  (go  to  (c)  above). 

4.  Calculate  the  source  term,  S,  Equation  (15),  then 
continue  iterating  velocity  potential  equation  (go 
to  1). 

The  number  of  inner  iterations  performed  on  the 
velocity  and  vector  potential  equations  are  arbitrary  and 
are  determined  by  experience.  Note  that  this  global 
iteration  strategy  has  not  been  optimized  in  terms  of 
convergence  rate  or  operation  count.  Its  primary 
purpose  is  to  provide  a  method  for  determining  the 
solution  to  the  coupled  system  of  equations. 

CONVECTION  OF  A  SINGLE  VORTEX 

The  convection  of  a  single  vortex  in  the  potential 
flowfield  about  an  axisymmetric  body  was  used  as  a  test 
problem  for  debugging  during  the  early  development  of 
the  computer  programs.  An  interesting  numerical  effect 
was  discovered  which  raises  some  questions  regarding 
the  accuracy  of  numerical  simulations  of  vortex  flows. 

'i 

A  conservation  equation  for  the  circulation  will 
first  be  derived  from  the  vorticity  transport  equation. 
Then,  an  equivalent  numerical  conservation  equation 
will  be  derived  from  the  difference  equations  used  to 
approximate  the  vorticity  transport  equation.  Finally, 
results  of  the  test  problem  will  be  given  to  illustrate  the 
degree  to  which  circulation  is  conserved  in  a  calculation. 
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If  the  vorticity  transport  equation,  Equation  (1),  is 
integrated  over  a  crossflow  plane,  one  obtains 


{  ||  di  +  J  $  •  (vx£nj  dA  -  0  (60) 


dA  -  dydz 


r*  V  -f  W  ^ 

Vvf  -  —  j  +  —  Jc 
00  00 


and  V  is  a  two-dimensional  gradient  operating  in  the  y 
and  z  directions.  Applying  the  divergence  theorem  and 
noting  that  A  *  f(x).  Equation  (60)  becomes: 


fc  J  fldA  + 

A 


8dl  -  0 


where  in  the  second  integral,  vn  8  is  the  vorticity  flux 
normal  to  1,  the  boundary  of  the  flowfield  in  the 
crossflow  plane.  Notice  that  the  first  integral  is  the 
definition  of  the  circulation,  r,  so  that  Equation  (61)  can 
be  expressed  as: 


F  =  Vfl  and  G  =  W8 

(The  last  term  on  the  right  hand  side  of  Equation  (63)  is 
the  approximate  factorization  error.) 

The  numerical  analog  of  the  circulation 
conservation  statement,  Equation  (62),  is  obtained  by 
summing  each  term  in  Equation  (63)  over  all  interior 
grid  points  in  the  respective  crossflow  plane.  If  the 
result  of  such  a  summation  involves  only  boundary 
terms,  the  numerical  method  is  said  to  be  conservative. 
This  scheme,  Equation  (63),  will  be  analyzed  below  to 
determine  if  it  is  conservative.  The  summation  of  each 
term  in  Equation  (63)  will  be  considered  separately, 
beginning  with  the  left  hand  side  term. 

Using  j  as  the  grid  point  index  for  the 
circumferential  direction,  rj,  and  k  as  the  index  for  the 
normal  direction,  $,  8i  + :  is  summed  over  all  interior 
points  in  the  i+1  plane  giving: 


j max-1  k max-1 


j max-1  k max-1 


i+1 / j / k 


Z  4- 


i+l/j/k 


vR0dl 


Equation  (62)  indicates  that  if  there  is  no  vorticity  flux  at 
the  boundaries  of  the  flowfield,  then  the  circulation  is 
conserved.  The  algorithm  used  to  solve  the  vorticity 
transport  equation  is  next  analyzed  to  determine  to  what 
extent  circulation  is  numerically  conserved. 

The  finite  difference  equations  which  determine 
the  convection  of  vorticity  over  one  step,  in  the 
streamwise  direction  are  represented  by  Equation  (55). 
Expanding  and  rearranging  Equation  (55)  gives: 


Mgi  +  Gi+i. 


#5  #6 

D.  +  {«  Vi  e.  -  6  V,  6  G  } 

i  4  (  i)  i  $  i  tj  i+1  $  i+lj 


Limits  on  the  summation  over  j  indicate  that  the  grid  is 
periodic  in  this  direction.  The  inverse  of  the  Jacobian,  J' 1 
represents  the  elemental  area  in  the  crossflow  plane 
associated  with  a  grid  point,  thus  J'1  =  AA  so  that 
Equation  (64)  can  be  rewritten: 


<«AA,.1+1(j(kC65) 


(The  limits  on  the  summation  in  Equation  (65)  and  all 
following  equations  in  this  section  are  implied  to  be  the 
same  as  in  Equation  (64),  i.e.  to  include  all  interior  grid 
points)  The  right  hand  side  of  Equation  (65)  is  a 
numerical  approximation  to  J8dA,  which  as  noted 
above,  is  equal  to  the  circulation,  f .  Consequently, 


Z  Z  « 


and  therefore  the  summation  of  + 1  over  all  interior 
grid  points* represents  the  circulation  contained  in  the 
i+l-th  crossflow  plane.  Similarly,  r.  is  obtained  by 
summing  8-  over  all  points  in  the  i-th  plane. 

Next  consider  the  summation  of  either  component 
of  term  #3  in  Equation  (63).  Making  use  of  the  difference 


where 
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formula  for  6  „  in  Equation  (53)  and  periodicity  in  the  j-  A  A  A  A  A  A 

direction,  it  can  be  shown  that:  -4ftj.t  +  6ftj-4ftjrt  +  ftj+2  (70a) 


z  z 

j  * 


5Fi, j,k  -  0 


(€7) 


The  difference  term,  S^F  is  therefore  conservative. 

For  6  ^ G .  in  term  #4  of  Equation  (63),  it  can  be 
shown  that: 


?  ¥  s^'k  “ 


z 

j 


{%i.i  + 


(68) 


^  ^Gi, j,kmax-l  +  Gi,j,kmax^  ) 

2  J 

Notice.that  the  two  terms  on  the  right  hand  side 
represent  approximations  to  boundary  fluxes,  located  at 
(i,j,k+l/2)  and  (i,j,kmax-l/2),  respectively.  Thus  in  the 
standard  nomenclature,  the  S^G  term  is  also 
conservative. 


It  can  also  be  shown  that  the  summation  of  the 
approximate  factorization  error,  term  tf 6  in  Equation 
(63),  is  identically  zero.  Using  Eqns.  (66)  -  (68),  the  result 
of  summing  Equation  (63)  over  all  interior  grid  points  is 

A  A 

r  -  r  +  M  v  I-  +  + 

i+1  l  2  ^  I  2 


+ 


*Gi,j,kmax-l  +  Gi,j,kmax) 
2 


and  " 


=  -4%-l  +  %'4«k+l  +k+2  (70b) 

i  A 

Note  that  for  applying  the  direction  term,  6^  ft, 
at  interior  grid  points  immediately  adjacent  to  a 
boundary,  a  formula  different  from  the  symmetric 
formula,  Equation  (70b),  must  be  used.  (Since  the  grid  is 
periodic  in  the  j-direction,  the  interior  point  formulation 
for  5*  ft,  Equation  (70a),  does  not  need  to  be  modified.) 
Various  boundary  point  formulas  were  tried  including  a 
second-order  difference  and  also  a  non-symmetric 
fourth-order  term.  However  with  these  formulations, 
summation  of  the  dissipation  term  results  in  a  non-zero 
linear  combination  of  near-boundary  ft  values  which 
represents  an  approximation  to  a  third  derivative  normal 
to  the  boundary.  This  boundary  term  can  become  quite 
large  when  there  are  significant  gradients  in  the  vortic- 
ity,  causing  the  terra  ZZD.in  the  numerical 
conservation  equation,  EquationT69),  to  act  as  a  source 
or  sink  of  circulation,  depending  on  the  term's  sign. 

To  remedy  this  problem,  the  following  boundary 
point  formulas  were  derived: 

sj a  k-2  -  A +  A  -  A +  «4  <7ia) 


3  A  A  A 

^  !k-kmax-l  ^kmax-3  ^kmax-2  + 

(71b) 

A  A 

3S*kmax-l  Scmax 


Using  Equations  (71)  with  Equation  (70),  results  in 


+  Z  Z  D± 

j  * 

It  is  thus  found  that  the  numerical  scheme.  Equation  (55), 
should  conserve  the  circulation,  provided  that  the 
dissipation  term,  Dj,  is  also  conservative.  This  term 
will  be  considered  next. 

The  dissipation  term  initially  used  in  the  program 
is  given  by  Equation  (56),  as  recommended  by  Steger 
and  Pulliam.6  However,  this  formulation  is  not 
conservative  because  the  inverse  Jacobian,  J’ 1 ,  appears 
outside  of  the  difference  operators.  This  .situation  is 
remedied  by  removing  the  inverse  Jacobian  and 
operating  directly  on  ft.  In  other  words: 


z  z  °ijk  -  o  <72> 

and  makes  the  dissipation  term  conservative.  By 
substituting  Equation  (72)  into  Equation  (69),  one 
obtains: 
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li+l  ‘i  +  2~  4-  l  2 


(73) 


*Gi, j,kmax-l  +  Gi, j,kmax* 


i“i,i+l 


This  final  equation  demonstrates  that  the  numerical 
approximation,  Equation  (55)  to  the  vorticity  transport 
equation  conserves  circulation  to  within  the  error  of  the 
boundary  flux  approximations. 

The  test  problem  used  to  demonstrate  the 
conservation  properties  of  the  numerical  method  will  be 
described  next.  An  axisymmetric  body  used  in  the 
simulation  consists  of  an  8-caliber  cylinder  having  a  3- 
caliber  tangent-ogive  cylindrical  nose.  A  base  having  the 
same  shape  as  the  nose  was  included  to  close  the  body. 
A  grid  was  constructed  about  this  body  having  41  x  21  x 
41  points  in  the  streamwise,  circumferential,  and  body- 
normal  directions,  respectively.  Views  of  the  grid  in  the 
symmetry  and  a  crossflow  plane  are  give  in  Figures  1(a) 
and  1(b),  respectively. 

Calculations  were  made  at  a  free  stream  Mach 
number,  Mtt  =  0.3,  and  angle  of  attack,  «  =  15  * .  The 
potential  flowfield  was  first  calculated  by  solving  the  full 
potential  equation.  A  single  vortex  was  then  placed  in 
this  flowfield  in  a  crossflow  plane  located  near  the 
beginning  of  the  straight  section  by  specifying  a  non¬ 
zero  value  of  vorticity  at  a  single  grid  point  in  the 
crossflow  plane.  Zero-normal  vorticity  flux  was 
specified  on  all  boundaries  in  all  crossflow  planes.  The 
vortex  was  convected  in  the  potential  flowfield  by 
solving  the  vorticity  transport  equation.  The  vortex  was 
advanced  twenty  streamwise  steps  downstream  to  a 
crossflow  plane  located  near  the  end  of  the  straight 
section.  Several  calculations  were  made  to  investigate 
the  effect  of  the  dissipation  on  the  conservation  of 
circulation. 

In  the  first  set  of  calculations,  the  vortex  was 
placed  at  a  location  approximately  one  body  diameter 
away  from  the  body  surface  in  the  initial,  upstream 
crossflow  plane.  This  put  the  vortex  twelve  grid  points, 
radially,  away  from  the  body,  indicated  by  the  symbol, 
®,  in  Figure  lb.  Two  calculations  were  made;  the  first 
using  the  non-conservative  dissipation  term.  Equation 
(56),  the  second  using  the  fully  conservative  dissipation 
term.  Equations  (70)  and  (71).  The  calculated  circulation 
in  each  crossflow  plane  was  evaluated  by  summing  U 
over  all  interior  points  in  the  plane,  according  to 
Equation  (65).  The  ratio  of  the  circulation  in  the  final 
and  initial  crossflow  planes  was  evaluated  to  measure 
the  degree  of  conservation  of  circulation.  Due  to  the 
specification  of  zero-normal  vorticity  flux  at  the 


boundaries,  the  initial  circulation  should  ideally  be 
preserved.  The  results  are  presented  in  Table  1. 

With  the  non-conservative  dissipation  term.  Case 
No.  1,  there  is  an  incremental  loss  of  circulation  as  the 
vortex  is  advanced  at  each  crossflow  plane,  resulting  in  a 
significant  net  loss  (14%)  ift  circulation  after  the  vortex 
has  been  advanced  to  the  end  of  the  straight  section. 
Using  the  conservative  dissipation  term,  the  loss  in 
circulation  is  significantly  less,  only  2%  of  the  initial 
value.  However,  no  loss  in  circulation  is  expected  since 
the  numerical  method  has  been  shown  above  to  be 
"conservative."  This  apparent  anomaly  will  be  discussed 
below. 

Another  series  of  calculations  were  made  with  the 
location  of  the  initial  vortex  much  closer  to  the  body,  in 
this  case  only  four  grid  points  radially  away  from  the 
body  surface,  indicated  by  the  symbol,  X,  in  Figure  lb. 
Both  methods  show  a  much  larger  loss  in  circulation  in 
this  instance.  Case  No.  2,  with  a  54%  loss  with  the  non¬ 
conservative,  and  a  22%  loss  with  the  conservative 
dissipation  term.  The  latter  loss  in  circulation  is 
unexpectedly  large  for  a  numerical  scheme  which  is 
supposed  to  conserve  circulation. 

The  source  of  the  loss  in  circulation  was  found  to 
be  the  boundary-normal  vorticity  flux  term,  G.  ^  /2  • 


Gl,j,l  +  Gl,j,2 
2 


(74) 


which  appears  in  Equation  (73).  Recall  that  zero-normal 
vorticity  flux  is  specified  on  the  crossflow  plane 
boundaries,  so  that  Gt  .  ^  =  0  at  the  body  surface, 
defined  by  k  =  l.  (Similarly  Gj  ^kmax  =  0  is  specified 
at  the  far  field  boundary,  k  =  'kmax.)  However,  since 
Gj  j  2  is  a  calculated  value,  it  can  attain  a  non-zero 
value,  particularly  when  the  vorticity  is  convected  close 
to  the  body,  and  serve  as  a  boundary  flux  of  vorticity  out 
of  the  flowfield.  This  explains  why  the  loss  in  vorticity  is 
much  less  when  the  vortex  is  located  away  from  the 
flowfield  boundaries. 

Although  the  numerical  method  was  shown  to  be 
conservative  in  the  standard  sense,  the  occurrence  of  a 
non-zero  boundary  flux  with  zero  flux  specified  is  an 
artifact  of  the  finite  difference  scheme.  The  error  which 
is  induced,  the  loss  in  circulation,  can  be  minimized, 
although  not  totally  eliminated,  by  reducing  the  mesh 
spacing  in  the  ^-direction  at  the  boundary  surface. 
Alternatively,  the  problem  can  be  completely  avoided  by 
using  a  finite  volume  formulation,  where  boundary 
fluxes  can  be  exactly  speeded.  In  any  case,  for  flows 
which  critically  depend  on  the  vorticity  contained  within 
the  flowfield,  it  is  of  utmost  importance  that  the 
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numerical  method  most  nearly  maintain  conservation  of 
circulation  to  ensure  a  high  degree  of  accuracy. 

RESULTS 

This  section  presents  computational  results  made 
with  the  computer  program,  TSBVRT,  (transonic, 
slender-body,  vortical  flows)  which  uses  the  algorithm 
described  above.  Comparisons  are  made  with 
experimental' data.  Results  are  first  described  for 
subsonic  flow  over  a  cylindrical  body,  followed  by 
results  for  subsonic  and  a  transonic  flow  over  a  body 
having  an  elliptical-shaped  cross  section. 

The  body  for  the  subsonic  case,  an  ogive-cylinder 
with  a  3-caliber  nose  and  8-caliber  body  is  the  same  as 
described  in  the  previous  section.  The  flow  conditions 
are  also  the  same:  free  stream  Mach  number,  =  0.3 
and  angle  of  attack,  a  =  15  \  This  particular  model  and 
flow  conditions  were  selected  to  correspond  to  a  set  of 
available  experimental  data  (Ref.  7). 

The  grid  used  for  this  case  is  the  41  x  21  x  41  grid 
that  also  was  described  above  and  is  illustrated  in  Figure  ■ 
1.  Flow  calculations  were  performed  in  the  following 
manner.  First,  since  TSBVRT  did  not  have  the  Stratford 
separation  criteria.  Equation  (37),  programmed,  the 
separation  line  was  determined  from  a  calculation  made 
with  the  VTXCLD  program  (Ref.  8)  for  the  same  body 
and  set  of  flow  conditions.  VTXCLD,  a  discrete  vortex 
method,  compares  well  with  this  experimental  data  and 
also  has  the  Stratford  separation  criteria  built-in. 

Using  the  VTXCLD-determined  separation  line 
location,  the  vortical  flow  field  was  calculated  with 
TSBVRT.  The  vorticity  flux  at  the  separation  location  in 
each  crossflow  plane  was  specified  at  a  single  grid  point 
according  to  Eqn.  (36).  The  boundary  layer  edge 
velocity,  Vc  was  taken  to  be  the  potential  velocity  on  the 
body  surface  at  the  separation  location  in  each  respective 
crossflow  plane.  The  vortex  reduction  factor,  X  =  0.45 
was  found  to  give  the  best  agreement  with  the 
experimental  data. 

Calculated  surface  pressures  and  corresponding 
velocity  fields  at  two  crossflow  stations  are  shown  in 
Figures  2  and  3,  respectively.  On  the  abscissa  for  the 
pressure  plots,  0  =  0*  corresponds  to  the  windward 
meridian  and  0  =  180*  to  the  leeward  meridian.  Also 
plotted  in  these  figures  are  the  corresponding  pressures 
for  both  experimental  and  VTXCLD  results.  As  the 
pressure  curves  indicate,  the  current  results  are  found  to 
agree  fairly  well  with  both  experimental  data  and  the 
VTXCLD  calculation.  Note  in  the  velocity  plots  the 
vortical  flowfields  located  on  the  leeward  side  of  the 
cylinder.  The  close  proximity  of  the  calculated  vortices 
to  the  body  is  responsible  for  the  correct  modification  of 
the  pressure  distribution  on  the  upper  surface,  0  * 
120*  —180*. 


The  agreement  between  the  present  results  and 
vortex  cloud  is  sometimes  not  as  good  as  would  be 
expected  from  two  mathematically  similar  formulations 
as  can  be  seen  in  Figure  2.  A  test  of  both  methods 
indicated  that  the  dominant  cause  of  the  discrepancy  is 
the  computation  of  the  second  order  term  uR  which  is 
due  to  the  tilting  of  the  vortex.  A  comparison  of  this 
term  at  x/D  =  7.5  computed  by  both  methods  is  shown 
in  Table  2.  It  may  be  seen  that  at  the  location  of  the 
greatest  error  in  surface  pressures,  $  ~ 120*,  the  present 
correction  is  much  less  than  that  computed  by  the  vortex 
cloud  method.  This  error  occurs  below  the  vortex 
feeding  sheet.  The  second  order  correction  is  found  by 
an  integral  of  the  vortical  velocity  wR  .  If  the  vortex 
sheet  is  diffused  by  the  calculation,  the  integral  from  the 
outer  to  inner  boundary  may  be  expected  to  be 
independent  of  the  diffuse  nature  of  the  sheet  if  the  total 
vorticity  along  the  line  of  integration  is  conserved  since 
it  does  not  depend  on  the  details  of  the  flow,  only  the 
end  points.  Since  the  present  method  is  conservative 
and  integrates  the  vorticity  to  obtain  the  correction  there 
should  be  no  significant  error  due  to  numerical  diffusion 
of  the  vortex.  Consequently,  the  cause  of  the  difference 
is  not  clear  at  this  stage. 

The  second  body  investigated  in  this  study  has  a 
2:1  elliptical  cross  section  with  a  parabolic  planform. 
This  body  has  been  experimentally  investigated  (Ref.  9) 
at  subsonic,  transonic  and  supersonic  free  stream  speeds. 
A  grid  was  constructed  about  the  body  with  41  x  81  x  23 
points  in  the  streamwise,  circumferential,  and  body- 
normal  directions,  respectively.  Views  of  the  grid  in  the 
planform  plane  and  also  in  a  crossflow  plane  are  given  in 
Figures  4(a)  and  4(b),  respectively.  Note  in  Figure  4(a) 
that  the  body  is  closed  with  a  simple  circular  arc  base. 
Note  also  that  all  dimensions  have  been  normalized  to 
2  A,  the  maximum  major  axis  of  the  elliptical  body.  The 
total  length  of  the  body,  not  including  the  base,  is 
L/2A  =  4.3.  ■ 

The  first  calculation  made  is  for  a  subsonic  flow 
with  free  stream  Mach  number,  Mw  =  0.4  and  angle  of 
attack,  a  =  16.3  * .  As  was  done  with  the  circular  body, 
VTXCLD  was  run  first  to  obtain  the  location  of  the 
separation  line.  However  for  the  elliptic  body,  it  is 
relatively  easy  to  infer  the  separation  line  from  the 
experimental  data.  The  incipient  separation  plane  can  be 
estimated  by  assessing  where  the  measured  surface 
pressures  begin  to  vary  from  the  pressures  obtained 
from  a  potential  flow  computation.  Downstream  of  this 
location,  the  flow  separates  on  the  shoulder  of  the  body. 
The  separation  point  was  treated  as  above,  except  that  in 
this  case,  the  vortex  reduction  factor,  X  =  0.60  gave  the 
best  agreement  with  the  experimental  data. 

Surface  pressures  computed  at  several  axial 
stations  are  plotted  against  the  experimental  data  and 
also  the  VTXCLD  results  in  Figure  5.  On  the  abscissa  in 
these  curves,  4>  =  0*  corresponds  to  the  leeward  meridian 
and  4>  =  180*  to  the  windward  meridian.  Agreement 
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with  the  data  for  both  TSBVRT  and  VTXCLD 
computations  is  not  as  good  for  this  body;  however, 
TSBVRT  appears  to  capture  the  correct  trend  of  the  data. 
The  cause  of  the  discrepancy  between  predictions  and 
data  is  not  known  at  this  time.  Grid  refinement, 
especially  in  the  body-normal  direction  may  improve  the 
accuracy  of  the  calculation.  To  illustrate  the  vortical  flow 
field  computed  with  TSBVRT,  the  velocity  field  is  shown 
in  Figure  6  for  the  crossflow  plane  at  X/2A  =  3.1. 

The  second  calculation  made  with  the  elliptic  body 
is  a  transonic  flow  with  free  stream  Mach  number,  Mw  = 
0.95  and  angle  of  attack,  <*  =  16.9\  The  separation  line 
was  the  same  as  in  the  previous  calculation,  except  here, 
X  =  0.90.  The  computed  surface  pressures  at  several  axial 
stations  are  plotted  along  with  experimental  values  in 
Figure  7.  The  predicted  results  have  about  the  same 
degree  of  agreement  with  the  data  as  in  the  previous 
case. 

For  this  flow,  the  free  stream  Mach  number  is  not 
large  enough  to  produce  crossflow  shocks;  however,  a 
shock  does  occur  on  the  upper  surface  of  the  base. 
Figure  8  plots  vorticity  contours  just  upstream  and 
downstream  of  this  shock.  This  latter  result 
demonstrates  the  ability  of  the  current  method  to 
compute  vortical  flowfields  with  shock  waves. 

Concluding  Remarks 

The  equations  of  Klopfer  and  Nixon2  have  been 
coded.  The  computer  code  is  a  considerably  developed 
version  of  a  potential  code,  TWING.  Results  for  bodies 
at  high  angle  of  attack  and  at  subsonic  and  transonic 
Mach  numbers  have  been  obtained  and  compared  to 
experiment,  and,  for  subsonic  flow,  to  a  discrete  vortex 
method.  The  equations  of  Klopfer  and  Nixon  assume 
that  only  one  component  of  vorticity  is  important  but  the 
present  study  shows  that  the  vorticity  component  due  to 
the  vortex  rotating  away  from  the  body  must  be 
included.  This,  in  fact,  is  introduced  as  an  empirical 
convection  in  the  discrete  vortex  method.  In  the  present 
work  this  is  shown  to  be  a  second  order  correction  to  the 
basic  theory. 

There  is  considerable  discrepancy  between  the 
results  of  the  present  theory  and  both  the  discrete  vortex 
method  and  experimental  data.  Since  both  prediction 
methods  solve  the  same  equations  but  with  different 
algorithms,  it  suggests  that  the  flow  may  be  very 
sensitive  to  errors. 

It  has  also  been  found  that  numerical  errors  can 
dissipate  an  isolated  vortex  by  a  considerable  amount,  a 
fact  that  suggests  further  work  in  algorithm 
development.  It  is  suggested  that  the  use  of  a  finite 
volume  method  would  remedy  this  problem. 
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Table  1 

Convection  of  a  Single  Vortex 
Check  on  Conservation  of  Circulation 


Case  No. 

ft  initial  location 

"  rfinal/rinitlaP 

Non-conservative 

dissipation 

Conservative 

dissipation 

1 

1  diameter  from  body 
surface 

.14 

.02 

2 

close  to  body 
surface 

.54 

.22 

Table  2 

Second  Order  Correction  Term  at  x/D  =  75 


8 

VTXCLD 

90 

-0.0018 

100 

0.0031 

110 

0.0484 

120 

0.1095 

130 

0.1042 

140 

0.0931 

150 

0.0864 

160 

0.0677 

170 

0.0938 

180 

0.0999 

Present  Results,  Equation  (46) 

0.0 

0.0035 

0.0211 

0.0438 

0.0798 

0.1201 

0.1205 

0.1039 

0.0629 

0.0000 


1  O 
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Figure  5.  -  Surface  Pressure  for  Elliptic-Parabolic 
Body,  Mw  =  0.4,  a  =  16.3. 


Figure  6.  -  Vortical  Flow  Field  at  x/2A  =  3.1 
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Figure  8.-  Vorticity  Distributions  for  Elliptic- 
Parabolic  Body. 
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THE  EFFECT  OF  VORTEX  ROTATION  ON  THE  FLOW 
AROUND  A  MISSILE  AT  VERY  HIGH  ANGLES  OF  ATTACK 

David  Nixon* 

Nielsen  Engineering  &  Research,  Inc. 

Mountain  View,  CA  94043-2287 

Introduction 

The  problem  of  predicting  vortical  flow  around  a  missile  can  be  solved  by  using1  the 
Navier-Stokes  equations,  which  are  very  expensive  to  use,  or  by  simpler  methods,  such  as 
that  by  Mendenhall  and  Perkins2  for  subsonic  flow.  In  the  latter  method,  vorticity  is 
introduced  into  the  flow  at  specified  locations  on  the  body.  Both  the  location  and  the 
amount  of  vorticity  are  determined  empirically.  A  transonic  version  of  the  theory  of 
Reference  2  is  given  by  Klopfer  and  Nixon3  with  a  more  complete  theory  given  by  Nixon  et. 
al.4  Klopfer  and  Nixon3  used  slender  body  theory  to  simplify  the  Euler  equations  to 
approximate  the  flow  around  missiles  from  five  three-dimensional  partial  differential 
equations  to  two  three-dimensional  equations  and  one  two-dimensional  equation.  The 
dominant  feature  of  this  formulation  is  that  only  one  component  of  vorticity,  the  streamwise 
component,  is  important.  In  Reference  2  and  in  Reference  4,  it  is  shown  that  the  inclusion  of 
a  second  vorticity  component  in  the  analysis  is  critical  to  the  success  of  the  simplified 
method.  This  second  component  reflects  the  tilting  of  the  separated  vortex  in  the  vertical 
direction.  In  view  of  the  importance  of  this  second  component  of  vorticity,  the  theory  of 
Klopfer  and  Nixon3  is  extended  in  this  .tote  to  include  the  appropriate  higher  order  terms 
in  order  to  gain  further  insight  into  the  vortical  flow  around  missiles. 

One  of  the  conclusions  of  the  extended  analysis  is  that  vorticity  in  the  vortex  sheet  is 
not  aligned  with  the  local  streamlines,  consequently  producing  entropy.  In  contrast,  the 

i 

rolled  up  vortex  is  aligned  with  the  local  flow  since  it  can  produce  no  load.  At  higher  Mach 


*  President 


numbers,  if  the  vorticity  is  introduced  at  a  non-normal  angle  to  the  body  surface,  this 
entropy  effect  can  lead  to  a  significant  total  pressure  loss  in  the  rolled-up  vortex. 

Analysis 


> 


•o 


The  Euler  equations  for  steady  compressible  flow  are 

(pU)x+(pV)y  +(pW)z=0  (1) 

(pU2  +  p)x  +  (pUV)y  +  (pUW)2  =  0  (2) 

(pUV)x  +  (pV2  +  p)y  +  (pVW)z  =  0  (3) 

(pUW)x  +  (pVW)y  +  (pW2  +  p)z  =  0  (4) 

[(e  +  p)U]x  +  Ke  +  p)V)y  +  [(e  +  p)W]z  =  0  (5) 

and 

p  =  (y  - 1)  {e  -  p(U2  +  V2  +  Wz)/2)  (6) 


where  p  is  the  density,  U,  V,  and  W  are  velocity  components  in  the  cartesian  coordinate 
system  (x,  y,  z),  e  is  the  internal  energy,  and  p  is  the  pressure. 

Manipulation  of  the  Euler  equations  and  the  use  of  Gibbs  relation  leads  to  Crocco's 
equation 

q  x  H  =  TVS  (7) 

where  S  is  entropy,  T  is  temperature,  q  is  the  velocity  vector  given  by 

q=7U+7V  +  kW  (8) 

and 

+  <9> 

The  vorticity  vector  Cl  defined  by 

n  =  Vxq  (10) 

Equation  (7)  can  be  differentiated  to  give 

Vx(qxO)  =  VTxVS  (11) 


or,  using  Equation  (7) 
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V  x  (q  x  fi)  =  VT  x  (q  x  n)/T 
In  component  form  Equation  (12)  can  be  written  as 

uo,x  +  (vn,)y  +(WQ,)a  =n,uy  +  apx  + 

[Ty(U02  -  vn,)  -  Tz(wn,  -un^l/T 

+  Vfi2y  +(wa,)2  =OjVx  +  H3V-,  + 

[Tx(ua,  -vnjj-T^vnj  -wa,)]/T 

(U^)x  +  (v^y  +wn32  =n,wx  +a,wy  + 

[Tx(Uft,  -Wn,)-Ty(Wi^  -  Vfij)]/T 
where  Q  =7Qj  +"[02  +  kl^ 

At  high  angles  of  attack  the  flow  separates  and  a  vortex  sheet  feeds  two  concentrated 
vortices  as  shown  in  Figure  1.  In  the  previous  work  by  Klopfer  and  Nixon3  the  body  was 
considered  slender  with  a  typical  cross  section/length  scale  of  e.  At  very  high  angles  of 
attack  the  flow  is  similar  to  that  around  an  equivalent  body  represented  by  the  body  itself 
and  the  recirculating  region  bounded  by  the  feeding  sheets.  Clearly,  there  are  two  scales  in 
the  problem,  namely  the  body  width/length  ratio,  given  by  e  and  the  characteristic 
height/length  ratio  of  the  recirculating  region,  given  by  6.  The  length  scales  £  and  6  are 
shown  in  Figure  1.  The  dominant  geometrical  effect  of  the  separation  is  represented  by  the 
length  scales  s  and  S. 

Now  construct  inner  coordinates  by  letting 

x  =  x,  y  =  y/5,  z  =  z/e  (16) 

The  velocity  will  consist  of  the  freestream  components  Uc  and  VQ  with  perturbations  the 
order  of  the  scales;  thus 

U  =  U0  +  6Pu1  +  ePu2,  V  =  V0  +  6  v,  W  =  ew  (17) 

where  if  Q0  is  the  resultant  freestream  velocity  and  a  is  the  angle  of  attack 


(12)‘ 


(13) 


(14) 


(15) 
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Uo  =Qo  c0S0t'Vo  =Q0  sin« 


(18) 
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The  index  p  is  such  that 


p>l 


For  a  constant  energy  flow 

h  +  j-CU2  +  V2  +  W2)  =  constant  (19) 

where  h  is  specific  enthalpy.  This  suggests  that  since 

h  =  CpT  (20) 

where  Cp  is  the  specific  heat  at  constant  pressure,  the  temperature  T  may  be  expanded  as  a 
series  like 


T  =  TQ  +  6Tj  +  e%  +  0(62)  (21) 

where  the  form  of  the  velocity  expansions  of  Equation  (17)  has  been  used.  It  is  assumed  in 
this  analysis  that 

6  >  >  e  >  >  62  (22) 

which  reflects  the  experimental  observation1  that  the  vortex  tilts  upwards  more  than  it 
spreads  outwards.  It  is  shown  later  that  in  order  to  obtain  the  dominant  effects  of  vortex 
rotation,  terms  of  the  order  of  62  must  be  retained.  Expansion  of  Equations  (13-15)  to  this 
order  gives. 


5U  a.  +  [v  +  6v)n  ]  +  6(wOJ  + 

o  l-  o  1  -  1  - 

x  y  z 


c2  63  > 

r  Vvo  +  -  vV  n! 

z  z 1 


(23) 


«ep- V-  +  ~  ^  Tfl+  5Tl-Uon3  +  6T2-Uon3 
3  z  €  o  2-  2  z  y 

z 
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50on2-+  <Vo  +  6v)n2?  +  5<VrfV--  + 

X  z  2 


(24) 


2  2  3 

—  v-  a,+  s2a.v-  -  s2v  n.T„  -4-  t„  a,v  -  4-  vt„  a.  -  62vt.  n, 
e  z  3  lx  ol2e2-3oe  2-3  1-3 

x  z  z  z 


su  n.  +  t(v  +5v)n,]  +  6wa,  -  (62t.  u  + 

O  3-  O  3  -  3-  2X  O 

x  y  y 

(25) 

2 

ST  V  +  6  T0  V)a,  =  Sew-  n,  +  enow 
x-  o  2-  3  x  1  2  - 

y  y  y 

By  examining  the  dominant  terms  in  Equations  (24)  and  (25)  it  can  be  implied  that 


0,-0  (S2^) 

n3  -  o  (8finx) 


(26) 


Retaining  terms  of  0(S2)  and  using  Equation  (26),  Equations  (23)  -  (25)  can  be  written  as 


U<Ax  +  WVy  +  W<Vz  +  W  =  0  (27) 

Vo^2y  +(TXU0  +  Tzw) 0,  ^  (28) 

03  =  0  (29) 

Note  that  from  Equation  (26) 

(Ojl  »  1^2*  » 

which  indicates  that  the  vortical  components  are  separable. 


The  terms  involving  the  temperature,  T,  are  due  to  entropy  production.  If  there  was  no 
entropy  then  these  terms  would  vanish.  In  the  earlier  work  of  Klopfer  and  Nixon3 ,  entropy 
effects  do  not  appear.  This  is  because  the  dominant  vorticity  component  and  the  dominant 
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velocity  vector  are  aligned.  In  the  present  case,  the  second  order  equations  do  not  have  this 
property.  Because  entropy  is  produced  in  the  vortex  feeding  sheet  there  will  be  a  total 
pressure  loss  in  the  vortex.  If  the  only  vorticity  component  at  separation  is  ,  then  entropy 
production  there  is  zero  only  if  the  separation  line  is  a  crossflow  stagnation  point.  The 
initial  entropy  production  is  determined  by  the  relative  angle  between  the  separation  line 
and  the  local  velocity  vector.  In  the  present  model  it  is  assumed  that  the  vortex  sheet  is 
tangential  to  the  body  surface,  as  sketched  in  Figure  1;  This  implies  that  the  windward  side 
of  the  separation  is  not  a  crossflow  stagnation  point,  an  assumption  which  is  compatible 
with  the  physical  description  of  a  separating  boundary  layer.  If  the  separation  line  is  a 
crossflow  stagnation  point,  then  the  crossflow  velocity  term  V  cannot  be  split  into  the  form 
given  in  Equation  (17)  because  VQ  and  5  v  could  become  comparable  in  magnitude  but  of 
opposite  sign.  In  such  an  event,  the  present  perturbation  analysis  is  not  valid.  Because  of 
the  closer  approximation  of  the  present  flow  representation  to  the  physical  problem,  it  is 
suggested  that  entropy  production  should  be  included  in  the  separation  model  and 
consequently  the  conclusions  noted  above  are  valid.  The  most  important  of  these  is  that 
given  in  Equation  (30)  which  indicates  that  the  vorticity  components  have  separable  scales 
of  magnitude.  It  should  be  noted  that  the  temperature  effects,  which  represent  entropy,  are 
proportional  to  and  hence  become  more  important  at  higher  speeds.  This  in  turn, 
implies  that  the  effect  of  the  flow  on  the  relative  angle  of  flow  separation  may  become  more 
important  at  higher  Mach  numbers  because  of  the  entropy  contribution. 

Consider  now  the  case  of  incompressible  flow  in  which  the  entropy  effects  are  not 
present;  the  equations  approximately  governing  the  vorticity  transport  can  then  be  found 
from  Equations  (27)-(29)  to  be  * 


+  (V0j)y  +  (WOj)2  =o 

(31) 

< 

> 

II 

$ 

(32) 
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The  logical  choice  of  the  index  p  in  Equation  (19)  is  two  since  from  Equation  (33) 

03  =  Vx  -  Uy  -  0  (34) 

or,  in  inner  varables  to  a  first  approximation 

5v  -  6Pu*  -  -  0  (35) 

x  1  y 

which  requires  that 

p  =  2  (36) 

Consequently,  Equation  (35)  becomes 

vx=U7  (37) 

If  Equation  (32)  is  written  in  inner  variables  then  it  becomes 

V0‘¥y=5Vl  <38> 

From  Equation  (31)  it  may  be  seen  that  if  the  retained  terms  are  of  the  same  order  of 
magnitude  and  if  VQ  is  0(1)  then,  to  a  first  approximation 

<Vy  =Q(ta i>  (39) 

which  implies  that  is  independent  of  y  to  a  first  approximation.  Hence,  with  the  help  of 
Equation  (37),  Equation  (38)  can  be  integrated  with  respect  to  y  to  give  to  a  first 
approximation 

Von2  =  62u101  +  F(x,  z)  +  0(83)  (40) 

where  F(x,z)  is  a  function  of  integration. 

At  the  separation  line,  y  =  ys  c  p ,  u1 ,  and  H1  are  specified;  hence  Equation  (40) 


becomes 


(41) 


n2  “  n2sep 


ni(°o  +  6  u’) 


li  _  j-  V°o  +  6  V  j 


+  0(63) 


sep 


where  the  subscript  "sep”  denotes  a  value  at  separation.  Reverting  to  the  usual  physical 
coordinates.  Equation  (41)  becomes 


ni° 

°2  "  n2sep  +  V- 
o 


V 


sep 


(42) 


Equation  (42)  is  an  algebraic  relation  between  and  rather  than  a  partial  differential 
equation  and  the  induced  velocity  could  be  computed  fairly  easily.  Although  the  velocity 
due  to  can  be  derived  from  Equation  (42)  a  simpler  means  that  uses  Equation  (30)  is 
derived  in  Reference  4.  However,  the  present  formulation  leads  to  more  insight.  It's 
possible  that  an  extended  analysis  would  produce  a  similarly  simple  relationship  for  f^. 

The  extra  component  of  vorticity,  ,  is  derived  by  a  stretching  of  as  implied  by  the 
differential  equation.  Equation  (14).  Because  of  the  relation  given  in  Equation  (40)  the  total 
magnitude  of  the  vorticity,  O  is  given  by 


1/2 


-  in,  I  +  O(o) 10, |  (43) 


However,  the  vortex  is  inclined  at  an  angle  0V  given  by 


0 


arctan 


which,  from  Equation  (42),  is  given  by 


(44) 
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(45) 


T7  ^^1  )  /  3 

6v  =  arctan  [  _  -  _  ytiy  +  O  /nj  +  0(S  ) 
o  o  J sep 


Since  the  local  flow  vector,  0f,  is  given  by 

0f  =  arctan  (V/U)  (46) 

it  may  be  seen  that  in  general  the  vorticity  in  the  vortex  sheet  is  not  aligned  with  the  flow 
vector.  Although  the  above  analysis  is  for  incompressible  flow,  this  result  is  consistent  with 
the  argument  regarding  entropy  production  given  earlier. 

At  the  end  of  the  vortex  sheet  the  rolled  up  vortex  cannot  support  a  load  and  therefore 
must  be  aligned  with  the  flow  vector,  that  is 

0V  =  0f  (47) 


and  hence  to  0(6) 


V 

U  jc 


0, 


sep 


/n2  + 


n~  seP]  + 

1  c 


0(63) 


(48) 


where  the  subscript  "C"  denotes  a  value  at  the  rolled  up  vortex.  If  the  angle  of  the  rolled  up 
vortex  is  denoted  by  0C  then  it  follows  from  Equation  (48)  that 


tan0 

c 


^  {tan©c)~1R 
o 


(49) 


where 


R 


-n 


2  sep 


n 


'lc 


(50) 
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From  Equation  (30)  s  e  p  *s  niuc^  ^e5S  ^an  '  and  .since  U  and  VQ  are  of  the  same  order 
of  magnitude,  R  is  positive. 

Equation  (49)  can  be  solved  to  give 

tan9c  =  j-{-R  ±  [R2  +  4VC/V0]1/2}  (51) 

which  gives  a  relation  for  the  angle  of  the  rolled  up  vortex.  Since  0Q  mus.  zero  if  no 
vorticity  is  shed,  as  is  the  case  when  a  =  0,  it  follows  that  the  positive  sign  in  Equation  (51) 
must  be  taken.  Thus 

tan9c  =  H-R  +  [R2  +  4VC/V0]1/2)  (52) 

For  a  long  body  the  strength  of  the  rolled  up  vortex  will  increase  as  more  body  shed 
vorticity  reaches  it.  From  experimental2  evidence  the  angle  of  the  rolled  up  vortex  becomes 
constant  and  hence 


V  V  V  V 

tan0  -  r—  «  tana  +  0(6  )  -  constant  (53) 

c  u  V  u  V  *“ 

c  o  c  o 

since  Uc  =  UQ  +  0(62 )  from  Equations  (17)  and  (36).  Thus  Vc  /  VQ  is  approximately 
constant  and  hence,  from  Equation  (52)  R  must  be  constant.  If  the  shedding  rate  is  constant, 
that  is  (02  VQ  )s  e  p  and  (UO!  )s  e  p  are  constant,  it  follows  that  Q1  c  is  constant.  Because  the 
vorticity  in  the  rolled  up  vortex  accumulates  along  the  vortex  because  of  the  steady 
shedding  rate,  it  then  follows  that  the  vortex  must  'expand  in  size  as  the  distance 
downstream  increases  in  order  to  keep  Oic  constant.  This  is  in  accordance  with 
experimental  data. 

Concluding  Remarks 

i 

An  analysis  for  vortical  flows  that  includes  two  components  of  vorticity  has  been 
derived.  The  analysis  indicates  that  vorticity  in  the  vortex  sheet  does  not  align  with  stream 
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lines  if  it  is  introduced  at  a  non-normal  angle  to  the  body.  Because  of  this  fact,  entropy  is 
produced  which  will  cause  a  loss  of  total  pressure  in  the  rolled  up  vortex.  The  analysis  also 
indicates  that  the  rolled  up  vortex  will  widen  as  it  progresses  down  the  body,  provided  the 
shed  vorticity  is  constant,  a  fact  known  from  other  sources,  but  a  corroborating  factor  in  the 
present  theory.  A  critical  result  is  that  the  three  vorticity  components  are  separable  in 
magnitude  for  a  slender  body  leading  to  a  less  restrictive  theory  than  in  Reference  3.  Thus  it 
is  possible  that  for  the  flow  around  slender  missiles  at  very  high  angles  of  attack  a  much 
simpler  set  of  equations  than  the  Euler  equations  may  model  the  flow.  This  simpler  set  of 
equations  would  consist  of  two  three-dimensional  partial  differential  equations,  one  two- 
dimensional  partial  differential  equation  and  two  algebraic  equations. 
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1.  Introduction 

In  recent  years  it  has  been  possible  to  predict.the 
unsteady  transonic  flow  around  a  wing,  especially  those 
typical  of  commercial  aircraft,  in  a  fairly  efficient 
manner.  Frequently,  the  computer  codes  that  are  used 
are  based  on  potential  theory  and  arr  onsiderably  faster 
than  a  corresponding  calculation  tsing  the  Euler 
equations.  These  methods  have  been  developed  for 
aircraft  and  are  not  really  applicable  to  missile  flows 
where  the  effects  of  vorticity  due  to  flow  separation  are 
important;  the  potential  approximation  cannot  predict 
the  effects  of  vorticity  other  than  by  representing  a 
vortex  wake  by  an  infinitely  thin  sheet  which  is  excluded 
from  the  computational  domain.  This  model  is 
complicated  and  may  not  be  a  viable  option  for  routine 
calculations  around  real  aircraft  or  missile  configurations 
because  the  geometry  of  the  vortex  sheet  can  get  quite 
complex. 

Missiles  generally  have  low  aspect  ratio  fins  and 
slender  bodies,  neither  of  which  are  features  of 
commercial  aircraft,  and  the  airflow  is  dominated  by 
vortex  effects;  this  is  in  contrast  to  the  flow  over  a 
commercial  aircraft  which  is  attached.  For  steady  missile 
flow  a  variety  of  prediction  methods  are  available, 
ranging  from  panel  methods  with  the  addition  of 
nonlinear  vortex  dynamics  to  the  Euler  equations  or 
Navier-Stokes  equations.  The  Navier  Stokes  equations 
will  model  flow  separation  and  track  the  resulting 
vorticity  but  are  only  as  accurate  as  the  turbulence 
model  used  in  the  calculation.  This  is  critical  in  defining 
the  separation  line.  The  prediction  methods  based  on 
the  Euler  equations  need  some  empirical  criteria  to 
predict  the  location  and  strength  of  the  separating  vortex 
sheet;  the  Euler  equations  will  then  model  the  transport 
of  this  vorticity  reasonably  accurately. 

A  detailed  discussion  of  the  methods  used  to 
predict  steady  flows  over  missiles  is  given  in  Reference 
1.  In  order  to  predict  unsteady  flow  over  a  missile  to 
adequate  accuracy,  the  method  must  be  able,  therefore. 
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to  treat  unsteady  separated  flow.  Furthermore,  if  it  is 
necessary  to  predict  the  steady  flow  at  transonic,  rather 
.than  supersonic,  speeds  the  method  must  be 
computationally  efficient  since  unlike  supersonic  flow 
problems  which  can  use  a  marching  procedure,  transonic 
flow  methods  require  numerous  sweeps  of  the 
computational  domain,  leading  to  large  computational 
times. 

It  would  be  ideal  to  use  the  Navier-Stokes 
equations  to  model  unsteady  transonic  flow  around 
missiles,  but  there  are  several  difficulties  with  this 
approach.  The  most  obvious  difficulty  is  the  computer 
time  required,  which  is  several  orders  of  magnitude 
greater  than  that  required  for  a  potential  calculation. 
Even  if  a  dedicated  supercomputer  is  available  for  such 
calculations,  the  computation  of  the  necessary  flow 
separation  may  be  inaccurate  because  of  the  inherent 
inaccuracy  in  many  present  turbulence  models.  The  next 
best  approach  is  to  use  the  time  dependent  Euler 
equations  which  require  some  empirical  factor  to  initiate 
separation  and  compute  the  shed  vorticity.  This 
empiricism  would  be  a  development  of  the  steady  flow 
criteria  described  in  Reference  1.  If  the  separation  line 
and  shed  vorticity  can  be  predicted,  then  the  Euler 
equations  would  be  a  viable  model.  However,  since 
missiles  are  slender  some  further  approximations  can  be 
made.  This  paper  is  concerned  with  developing  a 
method  of  predicting  the  unsteady  pressure  distribution 
on  missiles  at  transonic  and  supersonic  speeds.  The 
concepts  developed  in  this  work  are  extensions  of  an 
earlier  analysis2  for  steady  flow.  The  approach  is  to 
make  as  much  use  of  the  existing  technology  as  possible. 
The  final  goal  is  to  develop  a  computer  code  capable  of 
predicting  the  unsteady  transonic  flow  about  missiles  for 
use  in  aeroelastic  calculations  or  during  maneuver.  The 
present  paper  is  an  account  of  the  second  phase  of  the 
work,  namely  the  preliminary  development  of  an 
unsteady  flow  variant  of  the  method. 

The.approach  is  based  on  experience  gained  in 
steady  subsonic  and  supersonic  flow  predictions  for 
missiles,  in  particular,  the  fact  that  if  a  separation  line 
and  the  strength  of  the  vorticity  introduced  at  that  line 
can  be  estimated  empirically,  the  governing  equations 
(such  as  the  Euler  equations)  will  represent  the  transport 


of  this  vorticity  accurately.  Because  the  computational 
time  for  the  Euler  equations  for  an  unsteady  calculation 
is  considerably  greater  than  that  for  a  potential 
calculation,  a  simplified  model  is  used. 

The  basic  equations  are  an  extension  of  those 
derived  by  Klopfer  and  Nixon3  for  steady  flow  and  are  a 
subset  of  the  Euler  equations.  For  slender  bodies,  the 
five  Euler  equations  are  reduced  to  an  equation  similar 
to  the  three-dimensional  unsteady  potential  equation,  a 
three-dimensional  vorticity  transport  equation,  and  a 
two-dimensional  Poisson  equation.  Apart  from  the 
reduced  set  of  equations,  a  second  advantage  is  that  one 
of  the  equations  is  almost  identical  to  the  potential 
equation  for  which  there  are  well  tested  computer  codes. 
This  allows  the  development  of  a  prediction  method 
based  considerably  on  proved  technology.  A  derivation 
of  the  equations  is  given  in  the  Appendix. 
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2.  Basic  Equations 

From  the  Appendix,  the  basic  equations  are  as 
follows: 

Vorticity  Transport 

C\  +  +  (vf))y  +  (wfl)z  =  0  (1 ) 

Rotational  Crossflow 

Ayy  +  Aa  =  -Cl  (2) 

Mass  Conservation 

pt  +  (pu)x  +  (pv)y  +  (pw)z  =  0  (3) 

where 


The  boundary  conditions  on  the  body  are  zero 
flow  through  the  body.  There  is  zero  flux  of  vorticity 
through  the  body  except  at  the  specified  separation  line. 

3.  Computer  Code  and  Numerical  Methods 

Equation  (6)  is  identical  in  form  to  that  solved  by 
Batina  et  al4  in  the  CAP-TSD  code.  This  code  is  capable 
of  predicting  the  unsteady  transonic  and  supersonic  flow 
over  an  aircraft  with  stores.  Consequently,  CAP-TSD  is 
used  as  a  basis  for  the  present  method.  In  thi.  initial 
phase,  only  a  body  or  fuselage  is  considered;  in  the  near 
future,  fins  will  be  added  to  the  model. 
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The  algorithm  used  in  CAP-TSD  is  approximate 
factorization  with  time  stepping  as  an  option;  this 
algorithm  is  retained  in  solving  Equation  (7).  A 
modified  version  of  this  algorithm  is  used  to  solve  the 
vector  potential  equation.  Equation  (2).  The  vorticity 
transport  equation  is  solved  by  using  another  modified 
variant  of  the  approximate  factorization  algorithm. 


<t>  is  the  velocity  potential,  A  is  the  vector  potential,  p  is 
the  density,  is  the  freestream  Mach  number,  n  is  the 
vorticity  in  the  streamwise  direction,  and  u,  v,  and  w  are 
velocities  in  the  coordinate  directions  x,  y,  z  respectively. 

Using  Equations  (3),  (4),  and  (5)  and  making  the 
usual  transonic  small  disturbance  (TSD)  approximation 
leads  to  the  following  modified,  TSD  equation. 


3.1  Continuity  Equation: 


Because  CAP-TSD  uses  a  shearing  transformation 
described  by:  £  -  £(x,  y);  rj  =  y;  £  =  z;  r  » t.  Equation  (6) 
is  reformulated  as: 
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The  time-marching  procedure  is  performed  using  a 
second  order  accurate  implicit  scheme,  thus 
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Linearizing  about  0n+1  =  and  recasting  the 

problem  in  approximate  factorization  form  yields: 
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3.2  Vector  Potential  Equation: 


Equation  (2)  does  not  have  explicit  time- 
dependence.  It  is  solved  at  each  time  step  and  in  each 
individual  streamwise  plane  by  using  successive  over- 
relaxation.  For  the  continuity  and  vorticity  transport 
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equations,  if  the  Newton  iteration  option  is  used,  then  A# 
and  AH  approach  zero  at  each  time  step,  and  the  solution 
is  time-,  le.  For  the  vector  potential  equation,  the 
solution  it.  aiso  time-accurate,  since  at  each  iteration,  AA 
can  be  made  arbitrarily  close  to  zero. 


3.3  Vorticity  Transport  Equation: 

The  time-dependent  vorticity  transport  Equation 
(1)  can  be  rewritten  as 
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The  corresponding  finite  difference  equation  is 
formulated  similarly  to  Equation  (10)  and  yields: 
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Using  Newton's  method  (0?+1  =  (1*  +  An),  and  expressing 
the  problem  in  approximate  factorization  form  leads  to 
the  equation: 
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4.  Boundaiy  Conditions 

The  solution  to  Equations  (1),  (2),  and  (3)  requires 
that  properly  formulated  boundary  conditions  be 
specified  for  0,  A,  and  0.  For  the  farfield,  these 
boundary  conditions  are  based  on  the  fact  that  the 
disturbance  velocities  must  vanish  away  from  the  solid 
boundary.  At  the  body,  the  flow  is  tangential  to  the 
surface,  and  a  normal  vorticity  jet  is  introduced  in  the 
inviscid  flowfield  at  the  point  of  flow  separation. 

With  CAP-TSD's  cartesian  grid,  the  farfield 
boundary  conditions  can  be  formulated  as: 


V0  =  0 


A  =  A„ 


0  =  0 


(20) 


About  the  symmetry  plane,  V  and  O  are  antisymmetric  in 
y;  therefore  at  y  =  0: 


0  y  +  Az  =  0 

(21) 

O  =  0 


The  main  problem  with  the  boundary  conditions  is 
the  representation  of  the  fuselage.  In  CAP-TSD  the 


-4- 


fuselage  cross-section  is  assumed  to  be  elliptic,  and  the 
zero  flux  boundary  condition  on  the  body  surface  is 
transferred  to  a  bounding  cartesian  grid  (the 
"computational  fuselage")  by  using  ideas  from  slender 
body  theory.  For  angles  of  attack  typical  of  missile 
flight,  however,  the  existing  treatment  of  the  fuselage  is 
inadequate  and  an  alternative  scheme  was  developed. 

In  this  modified  scheme,  the  assumption  of  slender 
body  flow  is  still  used,  in  particular,  the  consequence 
that  only  the  crossflow  is  important.  The  boundary 
conditions  on  the  cartesian  computational  boundary  are 
found  from  the  "change  in  thickness"  effect,  which  is 
relatively  accurate,  and  by  using  the  analytic  solution  for 
the  crossflow  around  an  ellipse  to  give  the  value  of  the 
normal  velocity  on  the  computational  boundary.  The 
magnitude  of  the  crossflow  is  U..  sina,  where  U„  is  the 
freestream  velocity  and  a  is  the  instantaneous  angle  of 
attack. 


ar  -  £  <5=5-  -  (24> 

where  crQ  =  *(Z0),  and  <r0  represents  its  complex 
conjugate. 

In  Equations  (23)  and  (24),  FP  and  FR  designate  the 
complex  potentials  representing  the  effects  of  potential 
and  rotational  components  of  the  flow  respectively. 

For  any  point  Z  along  the  top  or  bottom  of  the 
computational  fuselage,  the  normal  velocity  is  then 
prescribed  according  to: 

wP  "  Ra<~3§Edz)  -  +  sin<“>  <25> 

and 

w*  -  Re<"iETz>  '  <26> 


The  fact  that,  in  the  absence  of  vorticity,  the 
"standard"  form  of  the  TSD  Equation  (i.e.,  when  the 
product  0x0y  is  negligible  in  Equation  (7)  (no  swept 
shocks))  reduces  to  Laplace's  equation  in  the  crossflow 
plane  suggests  the  use  of  an  analytic  solution.  For  an 
elliptical  body  with  semi-axes  a  and  b  in  the  y  and  z 
directions  respectively,  this  analytical  solution  is 
obtained  by  using  conformal  mapping  of  an  ellipse  onto 
a  flat  plate.  Such  a  transformation  is  given  by: 

^  2  +  (Z2  -  a2  +  b2) 1/2 

a  +  b 

(22) 


Z  +  (Z2  -  a2  +  b2) 1/2 
with:  Z  =  -z  +  iy. 

In  cr-space,  the  elliptical  boundary  and  the 
symmetry  plane  (y  *  0)  are  mapped  into  the  real  axis. 
This  mapping  is  used  to  transfer  boundary  conditions 
from  the  true  body  surface  to  the  computational 
fuselage.  At  each  streamwise  cross-section,  the  model 
which  is  used  to  represent  the  flow  between  these  two 
surfaces  is  that  associated  with  a  point  vortex  in  a 
crossflow.  In  c-space,  the  irrotational  crossflow  is 
represented  by  a  uniform  stream  of  complex  velocity: 

-  -  .in  (a)  -J2-+-E-  (23) 

In  addition,  a  point  vortex  is  used  to  model  the 
rotational  component  of  the  flowfield.  Its  strength  T  is 
calculated  numerically  by  integrating  the  vorticity  in  the 
computational  domain,  and  its  location  in  Z-space  is 
determined  as  the  vorticity  centroid  (ZQ).  The 
corresponding  induced  velocity  at  point  <r  is  then  given 
by: 


Along  the  vertical  sides  of  the  cartesian  computational 
box,  the  spanwise  component  of  velocity  is  similarly 
given  by: 

v*  -  <27) 
and 

v*  -  tt>  ■  A>  <28) 

Equations  (25)  and  (27)  are  used  directly  to  supply  CAP- 
TSD  with  the  required  Neumann  boundary  conditions 
on  0.  On  the  other  hand.  Equations  (26)  and  (27)  specify 
the  tangential  derivative  of  A  along  the  computational 
boundary.  Therefore,  A  must  first  be  integrated  along 
the  boundary  to  provide  the  proper  Dirichlet  condition 
for  Equation  (2). 

The  addition  of  sin(a)  to  02  in  (25)  simply  reflects 
the  fact  that  in  CAP-TSD  the  grid  is  aligned  with  the 
freestream.  Consequently,  the  slender  body  is  in  effect 
considered  to  be  inclined  at  an  angle  (-a)  within  the 
computational  fuselage,  and  W  represents  the 
downwash  effect,  as  seen  from  a  referential  aligned  with 
0„.  This  configuration,  however,  leads  to  considerable 
inaccuracies  in  resolving  the  vortex  close  to  the  body. 
More  importantly,  in  the  physical  case  where  vorticity  is 
introduced  at  the  top  or  upper  right-hand  side  of  the 
cartesian  boundary,  fl  can  only  be  transported  outside 
the  computational  flowfield,  towards  the  body,  and 
cannot  be  resolved  explicitly. 

This  problem  is  circumvented  by  placing  the  body 
parallel  to  the  grid  and  considering  the  freestream  flow 
to  be  inclined  at  an  angle  of  attack  a.  As  a  first 
approximation,  the  disturbance  potential  0  can  be  solved 
using  Equation  (11),  as  long  as  a  remains  small,  although 
one  of  the  goals  of  the  present  study  is  to  push  the 
theory  well  beyond  its  limits,  i.e.,  for  finite  a.  in  the 
vorticity  transport  equation.  Equation  (20),  0  is  therefore 
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replaced  by  (0  +  sin(a)  z)  in  order  to  effectively  allow  the 
transport  of  vorticity  away  from  the  computational 
fuselage. 


numerical  implementation,  the  incremental  normal 
vorticity  flux  released  between  x  and  x  +  Ax  is  derived 
from: 


The  transfer  of  boundary  conditions  from  the  body 
surface  to  the  cartesian  computational  boundary  using 
potential  flow  around  an  ellipse  can  be  justified  on  the 
basis  of  the  following  arguments.  First,  the  use  of  an 
incompressible  boundary  condition  is  warranted  because 
the  crossflow  Mach  number  Mc  =  M„sina  is  order  a,  and, 
therefore,  compressibility  effects  in  the  crossflow  plane 
are  of  order  cr  and  can  be  neglected,  since  the  angle  of 
attack  can  be  initially  assumed  to  be  itself  of  order  e  (Ref. 
6).  Second,  the  strictly  two-dimensional  analysis  must  be 
valid  for  L  »  D,  i.e.,  to  the  order  of  slender  body  theory. 
Finally,  if  t  is  a  characteristic  time  scale  for  the  transfer 
of  boundary  conditions  from  the  true  body  surface  to  the 
computational  fuselage  (length  scale  D),  then  an  upper 
bound  for  r  is  given  by  t  =  D/(C(1-Mc)).  In  an  unsteady 
flow  of  reduced  frequency  k  =  wL/U^,  the  characteristic 
unsteady  time  scale  is:  T  =  2ffL/(UJ0,  and  the  ratio  of 
these  time  scales  is  given  by: 


T 

T 


s 
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thus  justifying  the  use  of  an  effectively  quasi-steady 
boundary  condition,  valid  for  reduced  frequencies  up  to 
order  unity. 

A  flow  separation  condition  is  simulated  by  the 
injection  of  vorticity  into  the  inviscid  flowfield  at  the 
point  of  separation.  For  steady  flow,  a  modified  version 
of  the  Stratford  criterion7  is  used  to  determine  the 
separation  line: 

, ,  dc'  1/2,  .0.1 

Cp8dsJ  K  x  10  "  0.35  sina  <3°) 


where  Cp'  is  the  modified  pressure  coefficient,  and  s  is 
the  virtual  length  of  the  turbulent  boundary  layer,  as 
seen  in  the  crossflow  plane8. 

The  strength  of  the  vorticity  jet  is  derived  from  the 
observation  that  for  a  flat  plate  boundary  layer  of 
thickness  6,  the  streamwise  vorticity  flux  per  unit  span  is 
given  by: 

6  2 

r  D 

I  VOdy  -  (31) 

O 

where  Uc  is  the  velocity  at  the  edge  of  the  boundary 
layer.  At  the  point  of  separation,  it  is  assumed  that  a  net 
fraction  X  of  this  vorticity  flux  is  injected  into  the 
frcestream.  This  method  has  been  formally  shown2  to  be 
equivalent  at  subsonic  speeds  to  the  "vortex  cloud" 
method  used  by  Mendenhall  and  Perkins?  In  its 


x+Ax  a2 

V  fids  -  iw  2 A  Ax  (32) 

n  2  sp 

'sp"  denotes  a  value  at  separation, 
and  S|  and  s2  are  values  of  the  curvilinear  coordinate 
along  the  body  in  the  crossflow  plane,  placed  on  either 
side  of  the  separation  point.  The  empirical  "vortex 
reduction  factor"  X  determines  the  amount  of  vorticity 
shed  at  that  point  and  is  set  to  be  equal  to  0.6. 

The  transfer  of  the  vorticity  jet  boundary  condition 
to  the  computational  box  requires  further  analysis.  For 
consistency  with  the  model  used  to  transfer  velocity 
boundary  conditions,  the  Stratford  criterion  can  be 
implemented  by  using  a  transformation  which  inverts 
Equation  (22),  to  obtain  Cp  at  the  body  surface  (see  Sec. 
6).  Knowing  the  resulting  separation  point,  Z#,  one  can 
obtain  the  value  of  the  stream  function  /  =  Im(F((/)},  at 
g*  =  x(Z').  The  location  of  the  vorticity  jet  at  the 
computational  boundary  could  then  be  determined  as 
that  which  satisfies  f  =  y*.  At  the  present  time,  the 
computer  code  is  set  up  to  prescribe  any  distribution  of 
normal  vorticity  fluxes  along  the  computational  fuselage 
boundary.  The  results  presented  in  this  paper,  however, 
consider  only  the  model  case  of  a  concentrated  vorticity 
flux  placed  at  the  top  surface  and  designed  to 
approximate  the  above  scenario. 

5.  Second  Order  Correction 

Nixon  et  al2  have  shown  that  a  second  order 
extension  of  the  theory  is  required  in  order  to  account 
for  the  vortex  tilting  phenomenon.  The  tilting  away 
from  the  axial  direction  is  obtained  by  retaining  the 
rotational  component  of  streamwise  velocity,  UK.  This  is 
achieved  by  integrating  the  spanwise  vorticity  in  the 
normal  direction,  which  yields  to  leading  order: 


z 


The  inclusion  of  this  higher  order  term  was  shown 
to  give  more  accurate  results  and  laid  a  theoretical 
foundation  for  the  "empirical  correction"  used  by 
Mendenhall  and  Perkins4. 

6.  Pressure  Coefficient  Calculation 

The  pressure  coefficient  at  the  surface  of  the  body 
is  calculated  using  the  isentropic  relation 


x  s. 


where  the  subscript ' 
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where  q2  is  the  squared  modulus  of  the  total  velocity 
+  ^R/  including  the  second  order  correction  UR. 

Since  q2  and  30/3t  are  not  readily  available  at  the 
body  surface,  owing  to  the  use  of  the  cartesian 
computational  boundary,  it  is  necessary  to  use  an  inverse 
transformation  relating  30/dt  and  q2  at  the 
computational  fuselage  to  3^/3tls  and  q,.2  at  the  body 
surface. 

This  can  be  achieved  by  using  the  inverse 
conformal  mapping  transformation: 

„  _  (a  +  b)  ( <r  +  (<r2-  4)1/Z) 

2 - 4 - 

(35) 
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For  any  first  grid  point  away  from  the 
computational  boundary,  Z  is  projected  orthogonally 
onto  a  point  Zg  which  belongs  to  the  true  body  surface. 
From  Equations  (23)  and  (24),  the  crossflow  surface 
velocity  at  Zs  can  then  be  related  to  the  velocity 
calculated  at  Z  by  observing  that  in  <7-space: 


dFI  dFI 
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where  k  is  a  complex  transfer  function  given  by: 


2n  las-*0 


Instead  of  dF/dcrl  0,  the  actual  velocity  calculated  at  Z  is 
mapped  in  <r-space,  so  that: 


<-w-iv>§ 


+  *r  (ov<rs,  <r0,D 

Finally,  the  complex  velocity  at  the  surface  can  be 
evaluated  in  Z-space,  using  the  original  transformation 
expressed  in  Equation  (22): 


+  *:(<r,<7s,<r0,r) 


The  squared  modulus  is  then  calculated  by  assuming 
that  the  dominant  (streamwise)  component  of  velocity  is 
equal  to  that  calculated  at  Z,  hence: 


qs  ■  (1+W2  +  Vs  +  Ms  (40) 


Similarly,  the  complex  potentials  at  <rs  and  a  can 
easily  be  related  using  the  following  relation: 


l-^ln 
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Hence  30/3t  ls  is  not  equal  to  30/3t,  and  must  be 
modified  in  order  to  account,  in  particular,  for  the 
temporal  variations  of  r  and  which  can  be  the  source 
of  additional  phase  lags. 

7.  Results 

A  time-accurate  calculation  was  performed  for  a 
missile-shaped  body  of  11:1  aspect  ratio  and  2:1  elliptical 
cross-section.  A  planform  section  of  the  body  in  the 
symmetry  plane  is  shown  in  Figure  1.  The  body  was 
configured  at  a  mean  angle  of  attack  a  =  15°,  and  the 
freestream  Mach  number  was  set  to  =  0.3. 

The  mesh  size  was  45  x  25  x  54,  and  its  full  extent 
in  the  crossflow  plane  is  represented  in  Figure  2.  The 
location  of  the  computational  fuselage  is  indicated  by  the 
inner  boundary  of  the  grid,  within  which  the  elliptical 
cross-section  of  the  body  has  been  outlined.  A  vorticity 
flux  was  introduced  along  the  central  portion  of  the 
body  (0.26  s  X/D  s  0.78)  in  the  form  of  a  vertical 
vorticity  jet  (the  arrow  in  Figure  3).  Figure  3  also 
indicates  the  convention  used  for  the  definition  of  the 
roll  angle  S.  To  avoid  any  possible  confusion,  a  note 
should  be  made  that  throughout  the  remainder  of  this 
manuscript,  the  symbol  D  is  used  as  a  reference  length 
scale  and  refers  to  the  half-axis  of  the  elliptical  cross- 
section  in  the  spanwise  direction.  Similarly,  X,  Y,  and  Z 
are  used  to  refer  to  the  streamwise,  spanwise,  and 
vertical  coordinates  respectively. 

In  order  to  verify  the  accuracy  of  the  analytic 
"transfer  of  boundary  conditions"  derived  in  Section  4 
(Equations  (25)  through  (28)),  pressure  distributions 
were  computed  for  a  purely  irrotational  calculation. 
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Figure  4  shows  Cp  distributions  as  a  function  of  the  roll 
angle  0,  taken  at  the  mid-chord  of  the  body.  The 
pressure  distribution  for  the  2:1  ellipse  corresponds  to 
the  body  previously  described,  whereas  the  case  labeled 
'l:!'  corresponds  to  an  equivalent  body  with  a  circular 
cross-section.  In  the  latter  case,  the  Cp  distribution 
reveals  the  expected  cosine-shaped  profile,  with  a 
minimum  of  Cpmin  =  -0.29,  which  is  the  value  expected 
from  potential  flow,  since  the  cross-flow  at  0  =  n/2  is 
q*  =  (2tana)2  =  0.287. 

Converged  time-accurate  solutions  were  also 
computed  with  flow  separation  (i.e.,  with  the  vorticity  jet 
turned  on).  Vorticity  contours  in  the  crossflow  plane  are 
shown  at  successive  locations  along  the  body 
(Figures  5(a)  through  5(d)),  for  a  steady  calculation  at 
a  =  15°.  These  indicate  the  occurrence  of  vortex  roll-up 
and  intensification  of  the  vortex  strength  with 
downstream  distance.  Further  examination  of  the  cross- 
sectional  vorticity  contours  and  of  similar  contours  in  the 
x-z  plane  indicates  that  the  lowest  vorticity  levels  tend  to 
propagate  in  space  at  an  angle  of  approximately  15° 
from  the  body  (i.e.,  similarly  to  passive  advection).  The 
larger  vorticity  levels,  on  the  other  hand,  remain  close  to 
the  body  surface,  as  expected  from  the  induction  due  to 
the  image  vorticity.  This  behavior  is  consistent  with 
experimental  observatipn. 

The  decomposition  of  the  flowfield  into  rotational 
and  irrotational  velocity  fields  is  illustrated  in  Figure  6. 
Figure  6(a)  shows  the  potential  velocity  component  (i.e., 
+  sin(a)k,  where  k  is  a  unit  vector  in  the  positive  z- 
direction)  at  X/D  =  0.67.  The  rotational  component  Vx  A 
is  represented  in  Figure  6(b)  and  illustrates  the  presence 
of  a  formed  vortex.  The  total  velocity  field  (Figure  6(c)) 
is  obtained  by  the  superposition  of  rotational  and 
irrotational  flow  components.  The  corresponding  spatial 
evolution  of  the  flowfield  along  the  body  is  documented 
in  Figure  7,  illustrating  vortex  formation,  strengthening, 
and  lift-up  away  from  the  body  surface. 

A  comparison  of  pressure  coefficient  distributions 
with  and  without  the  introduction  of  vorticity  into  the 
flowfield  is  provided  in  Figure  8  for  various  downstream 
locations.  Upstream  of  the  location  where  separation 
first  occurs,  irrotational  and  vortical  Cp  distributions  are 
virtually  indistinguishable  from  one  another,  as  shown 
in  Figure  8(a)  for  X/D  =  0  23.  This  location  incidently 
corresponds  to  the  end  of  the  nose  section  of  the  body 
(see  Figure  1).  Therefore,  the  asymmetric  pressure 
distribution  is  a  consequence  of  the  forebody 
configuration.  At  X/D  =  0.34  (Figure  8(b)),  the  weak 
rotational  component  of  the  flow  produces  a  net 
reduction  of  surface  velocities  in  the  vicinity  of  the 
separation  point  and  tnerefore  an  increase  in  pressure  on 
the  lee-ward  side.  This  results  in  a  local  reduction  in 
cross-sectional  lift.  The  effect  becomes  more  pronounced 
as  the  strength  of  the  vortex  increases  (Figure  8(c)). 
However,  at  sufficiently  large  downstream  distances 
(Figures  8(d)  and  8(e)),  the  strength  of  the  vortex 
becomes  such  that  high-speed  reverse  flow  occurs  at  the 
top  surface,  resulting  in  low  pressures  and  increased  lift. 


The  full  distributions  of  cross-sectional  lifts  and  pitching 
moments  are  compared  to  a  purely  irrotational 
calculation  in  Figures  9(a)  and  9(b)  respectively. 
Although  incipient  separation  takes  place  at  X/D  =  0.26, 
the  enhanced  lift  due  to  flow  separation  is  only  apparent 
for  X/D  a  0.4,  in  accordance  with  the  above  mentioned 
threshold  in  vortex  strength.  The  total  lift  coefficient  for 
this  case  was  found  to  be:  CL  =  1.20.  For  reference,  the 
equivalent  lift  calculation  using  potential  theory  without 
the  inclusion  of  vorticity  effects  yields  CL  =  0.54. 
Similarly,  the  total  pitching  moment  taken  about 
Xrcf/D  =  0.45  is  CM  =  -0.07  (versus  -0.25  for  a  purely 
potential  calculation). 

Hence,  the  inclusion  of  the  rotational  flow  is  seen 
to  be  responsible  for  drastic  differences  in  the  flow- 
induced  loads  on  the  body. 

The  results  described  in  Figures  3  through  9 
correspond  to  a  steady  configuration  at  a  =  15°.  The 
following  results  examine  the  case  of  flutter  about  a 
mean  angle  of  attack.  The  body  motion  i*» 
approximated  to  leading  order  by  imposing  an 
oscillating  freestream  in  the  crossflow  plane  equal  to 
U„sin(a),  for  which  a  describes  a  simple  harmonic 
motion:  a(t)  =  a0  +  aTsin(kt),  where  aQ  =  15°  is  the  mean 
angle  of  attack,  =  3°  represents  the  motion  amplitude, 
k  =  2^fL/U^  is  the  reduced  angular  freqrency  based  on 
body  length,  and  t  is  non-dimensional  time.  For  most  of 
the  unsteady  results  presented  in  this  manuscript,  k  is 
equal  to  In,  corresponding  to  a  flutter  period  of  unity. 
For  M„  =  0.3  and,  say,  L  =  10  m,  the  dimensional 
frequency  is:  f  =s  10  Hz,  which  is  within  the  range  of 
interest  for  aeroelastic  computations  and/or  active 
control. 

The  vector  plots  of  Figure  10  illustrate  the  motion 
experienced  by  the  vortex,  as  an  entire  cycle  of  the 
oscillation  is  completed.  At  phase  kt  s=  n/2  (Figure 
10(a)),  the  angle  of  attack  is  maximum.  In  a  quasi-steady 
sense,  the  strength  of  the  vortex  would  also  be 
maximum.  This  is  not  the  case,  however,  and  the  vortex 
continues  to  grow  in  strength  during  the  'pitch-down7 
part  of  the  motion  (Figure  10(b)).  This  is  because  the 
excess  vorticity  accumulated  upstream  takes  a  finite  time 
to  convect  to  the  current  location  of  X/D  =  0.67.  As  the 
motion  proceeds  (Figures  10(c)  and  (d)),  the  vortex 
strength  diminishes.  A  strong  correlation  can  also  be 
noticed  between  the  strength  of  the  vortex  and  its 
spanwise  location,  as  expected  from  the  induced  velocity 
due  to  its  image. 

The  corresponding  pressure  coefficient 
distributions  on  the  body  surface  are  shown  in  Figure  11. 
These  indicate  the  sensitivity  of  the  vortex  strength  and 
position,  as. well  as  the  associated  dynamic  loads  to  the 
unsteady  motion.  In  particular,  the  largest  negative 
pressures  on  the  lee-ward  side  of  the  body  (6  1 120° ) 
actually  occur  around  kt  =  n,  i.e.,  at  the  half-point  of  the 
downstroke.  The  global  effects  of  unsteadiness  on  cross- 
sectional  lift  distributions  can  be  assessed  from  Figure 
12.  At  X/D  a  0.67,  the  lift  is  maximum  at  kt  =  n,  in 
accordance  with  the  results  of  Figure  10.  The  crossing  of 
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the  lines  also  indicates  that  the  value  of  kt  associated 
with,  e.g.,  maximum  lift  may  vary  with  downstream 
position.  Finally,  it  is  interesting  to  notice  that  for  kt  =  0 
and  kt  =  tt,  the  lift  distributions  are  almost  identical 
upstream  of  flow  separation  (i.e.,  X/D  <  0.26),  whereas 
they  differ  greatly  when  vortical  effects  are  present. 

The  total  oscillating  lift  and  pitching  moment  were 
calculated  by  integrating  cj  and  c^.  Figure  13  represents 
the  instantaneous  angle  of  attack  and  lift,  as  a  function  of 
non-dimensional  time,  for  k  =  2tt.  Substantial  phase  lags 
(approximately  jt/2  for  the  lift)  with  respect  to  the 
motion  are  observed.  Similarly,  the  pitching  moment 
was  calculated  to  lag  the  motion  by  approximately  5r/6. 
Figure  14  presents  the  lift  coefficient  time  history  for  1.6 
periods  of  the  oscillation,  at  the  same  reduced  frequency, 
for  the  transonic  case:  M„  =  0.9.  The  solution  was  started 
from  rest.  It  illustrates,  as  in  the  results  of  Figure  13,  the 
finite  physical  time  required  for  the  vorticity  to  build  up 
and  produce  the  enhanced  lift 


q2  =  U2  +  V2+W2 

where  p  is  the  density,  U,  V,  and  W  are  velocity 
components  in  the  cartesian  coordinate  system  (x,y,z),  h 
is  the  specific  enthalpy,  and  p  is  the  pressure. 

Manipulation  of  the  Euler  equations  and  the  use  of  Gibbs 
relation  leads  to  Crocco's  equation 

Vh°  +  +  q  x  0  -  T  V  S  (7) 

where  S  is  entropy,  T  is  temperature,  q  is  the  velocity 
vector  given  by 

q=TU+7V  +  kW  (8) 

and 


8.  Concluding  Remarks 

A  theory  was  developed  to  treat  flow  separation 
and  related  vortex  effects  in  unsteady  transonic  flow 
around  slender  bodies.  This  theory  involves  the 
simultaneous  solution  of  a  modified  TSD  equation,  a 
vector  potential  equation,  and  a  three-dimensional 
unsteady  vorticity  transport  equation. 

The  implementation  of  the  theory  was  performed 
using  a  modified  version  of  the  CAP-TSD5  computer 
code.  This  modified  versaon  yields  convergent  and 
time-accurate  solutions.  It  is  shown  for  the  first  time  that 
realistic  high  angle  of  attack  configurations  may  be 
calculated  using  CAP-TSD,  thus  showing  considerable 
potential  for  aeroelastic  computations  and  unsteady 
aerodynamics. 

Appendix 

Derivation  of  Equations 

Analysis: 

The  Euler  equations  for  a  steady  compressible  flow 
are 


P,  +  (pU)x  +  (pV)y  +  (pW)j  =  0 

(1) 

(pu), + (pU2 + p)„ + (pUV)y + (pUW)z = 0 

(2) 

(pv),  •  (pUV)x  +  (pV2  +  p)y  +  (pV W)z  =  0 

(3) 

(px)t  +  (pUW)x  +  (pVW)y  +  (pW2  +  p)z  =  0 

(4) 

|p(h+l/2q2)-p][  +  [pU(h+l/2q2))x  + 

|pV(h+l/2q2)]y  +  |pW(h+l/2q2)]z  =  0 

(5) 

and 


*  - r  rs +  ?rj +  *  rs  (9> 

h0  is  the  stagnation  enthalpy.  The  vorticity  vector  n  is 
defined  by 

fi  =  V  x  q  (10) 

Equation  (7)  can  be  differentiated  to  give 


V  |3- +  v  x  (q  x  fi)  -VTx  VS  (11) 


or,  using  Equation  (7) 

|rn  +  V(qxft)  =  VT(q  X n)/T  (12) 

It  may  be  shown  that  to  a  first  approximation.  Equation 
(12)  can  be  written  in  component  form  as 

nlt + unlx + (VOj)  y + (wrty2 = ryjy  +  OjUz  (13) 

+  (Uf^x  *  Vfl2y  +  (Wn^z  =  a1Vx  +  Cl^/Z  (14) 

% + (UOj)x  +  (Vn 3>y  +  wa*z = n1wx + ryVy  (15) 

where 

O  =  iflj  +  j  O2  +  kflj 

Assume  that  the  vorticity  is  produced  on  a  slender 
body  where  the  thickness  to  length  ratio  is  characterized 
by  the  small  parameter  e.  Thus,  the  dimensions  of  the 
body  in  the  y  and  z  directions  are  of  order  £.  In  order  to 
make  the  dimensions  of  the  body  equal,  the  following 
transformation's  used. 
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X  “  X 

y  -  y/c 
z  -  z/c 


(16) 


equations  then  become  (dropping  the  superscript  tilde  in 
the  following  for  convenience) 


V,2- Wly  =  0;  V^-W^-Q, 

(23a) 

U,y-Vix  =  0;  Ury-v„  =  n2  =  0 

(23b) 

U|z-Wlx  =  0;  Utz-Wrx  =  n3  =  0 

(23c) 

In  addition  it  is  assumed  that  the  velocity 
components  U,  V,  W  can  be  expanded  in  the  usual 
slender  body  expansion  to  give 


Iq  Equations  (23b)  and  (23c)  the  equations  for  the 
rotational  components  simply  duplicate  the  irrotational 
component,  which  suggests  that  a  velocity  potential,  tf, 
exists  such  that 


U  -  U  (1  +  c2u) 

oo 

V  -  £  U  v 

OO 

W  -  £U  w 


(17) 


i  he  temperature,  T,  is  also  expanded  as  a  series;  thus 

T  =  Tjl  +  £mT1)  (18) 

where  mat.  Using  Equations  (16),  (17),  and  (18)  it  can 
be  shown  that  a  first  approximation  to  Equations  (13), 
(14),  and  (15)  is 


nIt  +  013?  +  (yCl{)9  +  =  0  (19) 

^2t  +  °2x  +  vfly"  +  (wr^z  55  tfliV*  +  f^Vj  (20) 

%  +  +  (wo^j  =  cf^w*  +  (^Wj?  (21) 

If  at  some  boundary  the  vorticity  that  is  initiated 
has  a  vector  in  the  x  direction,  then  Equations  (19),  (20), 
and  (21)  show  that  to  a  first  approximation  fi2  and  %  are 
negligible  in  comparison  with  which  is  then  given  by 

nu + nix  +  frnpy  +  (wn^  =  o  (22) 

Thus,  in  the  slender  body  approximation  one 
component,  the  crossflow  vorticity,  is  dominant  to  a  first 
approximation,  and  this  vorticity  is  transported 
throughout  the  fluid  without  interchanging  with  the 
other  components.  The  neglected  terms  are  of  the  order 
£0^  In  order  to  solve  Equation  (22)  it  is  necessary  to 
specify  the  boundary  conditions.  These  boundary 
conditions  are  the  location  of  the  separation  line  and  the 
magnitude  of  the  shed  vorticity.  These  must  be  found 
from  empirical  relations  such  as  those  used  by 
Mendenhall  and  Perkins  (Ref.  1). 

Assume  that  the  velocity  field  is  composed  of  an 
irrotational  part,  denoted  by  the  subscript  i,  and  a 
rotational  part,  denoted  by  the  subscript  r.  Assume  also 
that  only  the  component  of  vorticity  is  significant; 
that  is,  terms  of  order  are  negligible.  The  vorticity 


Ui  =  l  +  *x 

V,  =  *y  (24) 

W,  =  *2 

and  that 

Ur  =  0;  Vw  =  0;  Wrx  =  0  (25) 


A  vector  potential,  A  is  defined  as 

qr  =  Vx  A  (26) 

where 

^=7Ur+"jVr+kWr  (27) 

and 

A  =  iAj  + JA2  +  kA3  (28) 


Substituting  Equations  (27)  and  (28)  into  Equation  (23) 
gives 


and 


Since 


A,  +  A. 
lyy  Izz 

a2-o 

A3  -  o 


- 


Vr  =  Alz;  Wr  =  -A 


ly 


(29) 


(30) 


,  A2  =  A3  -  0 

the  subscript  "1"  will  be  omitted  in  the  following 
discussion. 

The  equations  governing  the  transport  of  vorticity 
are  Equation  (22)  and  the  following  equations 
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Ayy  +  Azz  “  “nl 


References 


(31) 


U  *•  1  +  0 

x 


V 

w 


p  +  A 

y  z 


(32) 


where  0  is  the  velocity  potential. 

From  Equations  (25)  and  (32)  it  follows  that 


The  standard  transonic  potential  wing  theory  can 
be  deduced  from  Equations  (1)  and  (2)  with  the 
irrotational  assumption  and  the  isentropic  relation 

p/py  =  p./pl  (33) 

Equation  (2)  can  be  written,  using  Equation  (1) 

Ut  +  UUX  +  VUy  +  WUZ  =  -Ipx  (34) 

Using  Equations  (32)  and  (34),  Equation  (35)  becomes 


0  +  (1+0  )  0  +  (0  +  A  )  0 

rwt  rx  rxx  ry  z  Yxy 


(35) 


+  (02  Ay)  0xz  -  Jy-i)  2  Ip 


P  1  V-l 


M  x 


Since,  from  Equation  (33),  A  is  not  a  function  of  x. 
Equation  (36)  can  be  integrated  to  give 


"p-l 


(y  -  dm: 

1  +  - 5 -  [-2*t  -  2*a 


(36) 


-  02  -  02  -  02  -  2A  0  +  2A  0 

jC  y  *  2  ^  y  y  z* 


J)'“ 


Equation  (37)  is  an  equation  for  p  in  terms  of  0,  A.  The 
set  of  equations.  Equation  (1),  the  irrotational  equations, 
and  Equation  (37)  with 

A2  =  Ay  =  0  (37) 

are  the  equations  solved  by  the  traditional  potential 
method.  In  order  to  solve  for  a  flow  with  vorticity,  two 
additional  equations,  namely  Equations  (22)  and  (31), 
must  be  solved.  Equation  (22)  gives  the  vorticity 
transport  and  Equation  (31)  the  rotational  velocity 
induced  by  the  vorticity. 
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Cross-sectional  View  of  the  Cartesian  Grid. 


Figure  4. 


Calculated  Irrotational  Pressure  Coefficient 
Distributions  for  1:1  and  2:1  ellipses. 


Cross-sectional  Vorticity  Contours  at  Various 
Stream  wise  Locations  (M^O.3,  cr=15°). 
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Figure  6.  Illustration  of  the  Helmholz  Decomposition 
Procedure  at  X/D=0.67  (M^O.3,  a=15°). 
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Comparison  of  Potential  and  Vortical 
Pressure  Coefficient  Distributions  at  Various 
Downstream  Locations  (M^O.3,  a=15°). 


9.  Comparison  of  Potential  and  Vortical  Cross- 
sectional  Lift  and  Pitching  Moment 
Coefficient  Profiles  (M^O.3,  a=15°). 
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APPLICATION  OF  INDICIAL  THEORY  TO  THE 
PREDICTION  OF  UNSTEADY  SEPARATION 
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1.  Introduction 

During  flutter  or  unsteady  maneuver  of  slender  bodies,  such  as  missiles  at  angle  of 
attack,  important  phase  delays  can  be  observed  between  the  unsteady  loads  acting  on  the 
body  and  the  body  motion  itself.  These  phase  delays  result  from  dynamic  effects,  and  in 
particular  the  fact  that  the  position  and  strength  of  vortices  can  lag  the  motion  of  the  body 
significantly.  One  of  the  principal  causes  for  this  effect  is  that  the  position  and  strength  of  a 
vortex  at  any  streamwise  location  along  the  body  are  the  integral  consequence  of  vorticity 
transport  dynamics  on  the  one  hand,  and  the  unsteady  rate  at  which  vorticity  is  fed  into  the 
vortex  sheet  on  the  other.  Therefore,  the  "history  effect"  may  be  thought  of  as  a  delay  due  to 
convection  and  to  the  separation  process  itself.  It  is  to  the  latter  that  this  manuscript  is 
devoted. 

The  problem  can  be  stated  as  follows:  given  a  small  change  in  instantaneous  "outer 
flow  conditions"  at  a  given  downstream  location  (these  conditions  already  include  the 
various  convective  phase  delays),  what  are  the  additional  delays,  if  any,  that  are  associated 
with  the  separation  of  the  boundary  layer? 

If  Ref  1,  the  essence  of  boundary  layer  separation  at  high  Reynolds  numbers  is 
characterized  by  a  normal  vorticity  jet.  With  this  representation  in  mind,  the  goal  of  this 
paper  is  to  establish  whether  the  unsteady  characteristics  of  boundary  layer  separation,  i.e., 
the  motion  of  the  separation  point  and  the  unsteady  vorticity  flux  shed  at  that  point,  can  be 
predicted. 

The  approach  is  to  analyze  the  motion  of  the  separation  point  as  well  as  the  vorticity 
flux  in  a  given  cross-flow  plane,  i.e.  on  a  circular  body,  using  two-dimensional  Navier- 
Stokes  calculations.  The  final  goal  is  to  incorporate  the  finding  of  this  study  into  a  modified 
version  of  the  Stratford  criterion,2  relating  the  separation  roll  angle  to  the  pressure  gradient 
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history  at  the  surface.  An  innovative  avenue  of  research  has  been  undertaken,  in  which  an 
attempt  is  made  to  predict  the  characteristics  of  unsteady  separation  based  on  the 
knowledge  of  the  impulse  response  of  the  separated  flowfield. 

2.  Indicia!  Theory 

*  « 

Assume  that  a  quantifiable  entity  /i(t)  characterizing  some  aspect  of  flow  separation  can 
be  defined  (in  our  case  //  might  typically  represent  the  separation  roll  angle  9S  or  the  normal 
vorticity  flux  associated  with  boundary  layer  separation,  but  is  not  restricted  to  these 
quantities).  Consider  a  hypothetical  change  in  p(t)  due  to  an  infinitesimal  change  in  some 
parameter  c(t)  at  some  time  t.  The  parameter  c  represents  a  change  in  outer  flow 
conditions  which  affects  separation,  e.g.  e  may  represent  a  small  change  in  angle  of  attack, 
or  a  small  change  in  crossflow  velocity.  Following  the  derivations  of  Reference  3,  if  p  varies 
continuously  with  c,  then  the  change  in  /i(t)  may  be  expanded  according  to 

5//(t)  -  /(e(t,T)  ^iili  6r  +  0 (6t ) 2  U) 

where  t  is  time,  5/j(t)  is  the  change  in  p  resulting  from  a  change  in  the  forcing  parameter  e  at 
some  time  t,  and  p^( t,r)  designates  the  rate  of  change  with  respect  to  e  of  pit)  at  time  r. 

Neglecting  the  higher  order  terms  in  Equation  (1),  the  total  integrated  effect  of  such 
steps  in  e  from  t  =  0  until  time  t  can  be  expressed  as: 

t 

A/i(t)  -  //£(t,r)e(0)  +  j  /<e(t,T)  dr  (2) 

o 

If,  in  addition  it  can  be  assumed  that  the  behavior  of  p  with  e  is  linear,  then  pe(t,r)  can  be 
represented  by  its  value  at  r  =  0  provided  that  t  is  taken  relative  to  x.  Thus 

//e(t,x)  a  fJe(t  -  T,0)  a  fl£(t  -  T)  (3) 

where  the  functional  form  of  p£  has  been  contracted  for  notational  convenience. 

Substituting  Equation  (3)  in  Equation  (2)  yields,  after  a  simple  change  of  variable 

t 

AMt)  “  He(t)c(0)  -  |  ne( T)  d£  dT  (4) 

o 

Equation  (4)  states  that,  under  the  linearity  assumptions  afore  mentioned,  one  can 
analytically  predict  the  change  A/i(t)  due  to  an  arbitrary  change  in  the  forcing  parameter  e, 
provided  that  the  indicinl  response,}! €,  is  known.  Equation  (4)  illustrates  the  power  of  the 
indicial  method,  since  the  response  A/i  may  be  calculated  from  the  knowledge  of  the 


-  ? . 


excitation  parameter  £,  and  that  of  the  indicial  response  /ig,  which  requires  only  a  single 
experimental  or  numerical  determination. 

Although  indicial  theory  has  traditionally  been  used  in  unsteady  aerodynamics  (see 
e,g.  Tobak/1  1954,  Jenkins,5  1988)  the  work  presented  in  this  manuscript  is,  to  the  authors' 
knowledge,  the  first  attempt  at  applying  this  concept  to  the  problem  of  flow  separation. 

The  approach  is  to  impose  a  step  change  in  Mach  number  (the  case  of  Reynolds 
number  changes  is  deferred  to  the  Appendix)  to  the  separated  flowfield  about  a  circular 
cylinder,  and  record  the  indicial  responses  of  separation  roll  angle,  shed  vorticity  flux,  and 
drag.  Once  these  indicial  responses  are  recorded,  the  validity  of  the  indicial  method  is 
tested  by  evaluating  the  accuracy  of  the  prediction  given  by  Equation  (4)  against  the  actual 
(numerically  computed)  change  A/i(t)  due  to  some  arbitrary  change  in  £.  Details  of  this 
procedure  are  given  in  the  Results  section. 

3.  Computer  Code  and  Boundary  Conditions 

The  two-dimensional  Navier-Stokes  calculations  were  implemented  on  a  245  x  79  C- 
mesh  grid.  The  grid  was  a  four-fold  spatial  refinement  of  that  used  by  Rodman,6  with  a 
symmetry  boundary  condition  on  the  centerline.  This  symmetry  condition  ('inviscid  splitter 
plate'  condition)  at  y  =  0  was  imposed  upstream  and  downstream  of  the  circular  cylinder  so 
as  to  prevent  vortex  shedding  beyond  the  critical  Reynolds  number.  This  was  a  necessary 
precaution  since  this  regime  is  not  observed  in  the  three-dimensional  flow  (i.e.  slender  body 
at  moderate  angle  of  attack)  of  which  the  two-dimensional  cylinder  constitutes  a  simplified 
crossflow  representation.  The  C-grid  had  a  total  streamwise  extent  of  25  diameters,  and 
extended  5  diameters  in  the  normal  direction,  away  from  the  centerline. 

A  detail  of  the  grid  is  shown  in  Fig.  1.  In  this  Figure,  the  flow  direction  is  from  left  to 
right,  and  the  sign  convention  for  the  roll  angle,  0,  is  taken  to  be  in  the  clockwise  direction, 
for  consistency  with  references  1  and  7.  (i.e.  0  =  0°  at  the  windward  stagnation  point).  The 
computer  code  that  was  used  for  all  of  the  calculations  presented  in  this  paper  was 
ARC2D.8  This  code  uses  an  approximate  factorization  finite  difference  scheme.  The  time 
integration  was  performed  with  a  3-point  second  order  accurate  implicit  method.  At  the 
outer  boundaries,  a  characteristic-like  boundary  condition  procedure  using  locally  one¬ 
dimensional  Riemann  invariants  (see  Pulliam,  1984)  is  implemented.  At  the  cylinder 
surface,  the  no-slip  boundary  condition  is  applied. 

In  order  to  simulate  instantaneous  changes  in  angle  of  attack  of  the  three-dimensional 
slender  body,  the  two-dimensional  cylinder  is  accelerated  through  the  flow  by  moving  the 
computational  grid.  The  ARC2D  code  possesses  time-metric  terms  which  allowed  the 
computational  grid  to  be  translated  with  relative  ease.  Using  a  minor  amount  of  post¬ 
processing  it  is  then  possible  to  examine  the  solution  flowfield  in  a  frame  of  reference 
attached  to  the  body,  therefore  simulating  a  change  in  freestream  velocity. 

For  a  viscous  compressible  flow,  conservation  of  momentum  may  be  expressed  as: 
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where  the  velocity  components  u.  are  non-dimensionalized  by  the  freestream  velocity  UTO. 
The  coordinates  Xj  are  non-dimensionalized  by  the  cylinder  diameter  D,  and  p,  p,  and  p  are 
non-dimensionalized  by  their  respective  freestream  values.  In  Equation  (5),  use  is  made  of 
the  perfect  gas  assumption 


P 

oo 


P  U2 

oo  oo 


(6) 


and  the  Reynolds  number  Re  is  based  on  the  freestream  velocity:  Re  =  PoPJO/p^.  Hence, 
with  this  formulation,  a  change  AU^  in  freestream  velocity  results  in  the  simultaneous 
change  of  both  Mach  number  (AMOT  =  AU^/a^)  and  Reynolds  numbers 

(ARe=p„AUJD/fO. 

It  is  possible  to  vary  only  one  parameter  if,  instead,  the  velocity  components  are  non- 
dimensionalized  by  the  freestream  speed  of  sound  a^.  Substituting  u.'  =  M^Uj  and 
t'  =  t/M^  in  Equation  (5)  and  using  continuity  yields  (dropping  the  prime  superscript  for 
notational  convenience). 
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is  the  Reynolds  number  based  on  freestream  sound  speed. 

Equation  (7)  is  precisely  the  form  of  the  momentum  equations  being  solved  by  ARC2D. 
Therefore,  for  a  fixed  cylinder  diameter  D  and  given  fluid  properties,  the  Reynolds  number 
Rea  remains  constant.  Hence  any  change  of  the  freestream  velocity  is  obtained  either  by 
changing  the  Mach  number  (which  then  affects  the  solution  through  the  farfield  boundary 
conditions),  or  by  translating  the  computational  grid  and  imposing  that  the  velocity  at  the 
cylinder  surface  equals  that  of  the  grid. 

Since  the  non-reflective  characteristic  outer  boundary  conditions  use  locally  one¬ 
dimensional  approximations,  it  was  determined  that  the  most  accurate  way  of 
implementing  a  change  in  freestream  velocity  is  to  set  the  computational  grid  in  motion. 
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An  example  of  this  procedure  is  shown  in  Figure  2.  A  steady  state  solution  with 
crossflow  Mach  number  Mc  =  0.25  is  first  computed.  At  t  =  0,  the  cylinder  is  impulsively 
translated  either  with  or  against  the  flow  (top  and  bottom  graphs  respectively),  at  a  speed 
equal  to  some  percentage  of  the  freestream  velocity.  Figure  2  illustrates  typical  velocity 
profiles  in  a  reference  frame  attached  to  the  cylinder,  at  a  time  tUc/D  =  1.0  after  the 
initiation  of  the  grid  motion.  The  top  graph  corresponds  to  an  effective  Mach  number 
change  AMC  =  -0.025,  while  the  bottom  graph  (cylinder  moving  against  the  flow) 
corresponds  to  AMC  =  +0.025. 

4.  Results 

In  order  to  simulate  the  steady  symmetric  leeward  vortex  system  present  in  slender 
body  flow  at  moderate  angle  of  attack,  the  asymmetric  vortex  shedding  phenomenon  was 
inhibited  by  imposing  a  symmetry  boundary  condition  at  y  =  0.  At  the  high  Reynolds 
numbers  typical  of  missile  flow,  the  separating  boundary  layer  takes  the  form  of  a  thin 
vortex  sheet.  The  separation  process  can  therefore  be  simulated7  using  a  concentrated 
normal  vorticity  jet  placed  at  the  body  surface.  In  a  two-dimensional  simulation  of  the 
crossflow,  however,  the  Reynolds  number  based  on  diameter  must  remain  low  (i.e.  of  the 
order  of  100  or  less)  in  order  to  best  approximate  the  topology  of  the  flowfield.  This 
constraint  originates,  from  the  fact  that  in  two  dimensions  vortex  shedding  occurs  beyond  a 
critical  Reynolds  number  Re  =  40.  If  shedding  is  suppressed  by  means  of  a  splitter  plate,  the 
flow  exhibits  elongated  regions  of  separated  flow  which  do  not  adequately  represent  the 
leeside  vortices  of  an  inclined  cylinder.  For  this  reason,  the  two-dimensional  Navier-Stokes 
calculations  presented  here  are  restricted  to  the  Reynolds  number  range  20  <  Re  <  90. 

At  these  relatively  low  Reynolds  numbers,  the  separating  boundary  layer  does  not  take 
the  form  of  a  vortex  sheet,  but  rather  that  of  a  thick  vorticity  layer  separating  from  the 
surface,  as  seen  from  the  vorticity  contours  of  Figure  3.  A  close  examination  of  these  typical 
vorticity  contours  reveals  that  the  vorticity  jet  representation  may  not  be  adequate  at  the 
low  Reynolds  numbers  which  are  needed  to  match  the  range  of  observed  separation  angles 
in  the  three-dimensional  configuration. 

These  observations  suggest  that  a  criterion  which  should  emulate  the  asymptotic  (high 
Reynolds  number)  vorticity  jet  representation  must  be  based  on  the  tracking  of  the  line  of 
vorticity  maximum.  Such  a  criterion  should  also  prove  to  be  superior,  in  principal,  to  the 
detection  of  flow  reversal  at  the  surface  since  the  latter  is  known  to  be  inadequate  in 
unsteady  flow  as  an  indicator  of  flow  separation  location. 

A  typical  plot  of  the  radial  location  of  vorticity  maximum  as  a  function  of  roll  angle  is 
shown  in  Figure  4.  It  can  be  established  from  this  plot  and  an  analysis  of  the  corresponding 
steady  flowfield  that  as  long  as  the  flow  remains  attached,  the  maximum  of  vorticity  lies 
along  the  surface  (i.e.  r/D  =  0.5  in  Figure  4).  At  a  critical  roll  angle  @s  (separation  roll  angle), 
the  vorticity  maximum  departs  from  the  surface.  In  the  asymptotic  high  Reynolds  number 
case,  the  locus  of  vorticity  maxima  presumably  coincides  with  the  location  of  the  separating 
vortex  sheet.  Based  on  these  remarks,  the  separation  roll  angle  was  defined  as  follows: 
according  to  the  previously  described  detection  criterion,  the  maximum  of  vorticity  is  found 
at  the  surface  when  the  boundary  layer  is  attached.  Hence  in  that  region,  the  normal 
gradient  80/  8 n  (where  n  designates  the  outward  normal  coordinate)  is  negative.  On  the 
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other  hand,  when  the  flow  is  separated,  the  vorticity  maximum  occurs  away  from  the 
surface  (r/D  >  0.5)  and  the  normal  vorticity  gradient  30/ 3n  becomes  positive. 


Therefore,  the  zero-crossing  of  the  normal  vorticity  gradient  can  be  found  by  linear 
interpolation  of  30/ dn  with  respect  to  0,  and  subsequently  used  as  a  continuous  indicator 
of  the  separation  roll  angle.  The  linear  interpolation  scheme  also  allows  to  partially 
overcome  an  inherent  lack  of  spatial  resolution,  since  changes  in  the  separation  location  0S 
are  often  quite  small  and  typically  of  the  order  of  the  grid  resolution  itself  (see  Appendix  A). 
As  mentioned  in  Section  2,  6S  thus  defined  represents  one  of  the  quantifiable  entities  /i(t) 
characterizing  flow  separation.  The  second  quantity  of  interest  is  the  vorticity  flux  at  the 
point  of  separation.  The  Stratford  criterion2  assumes  that  a  fraction  X  of  the  tangential 
vorticity  flux  across  the  attached  boundary  layer  gets  injected  into  the  freestream  at  the 
point  of  separation.  It  is  justified,  therefore,  to  consider  the  time-dependent  behavior  of  the 
tangential  vorticity  flux  across  a  fixed  contour  which  cuts  the  boundary  layer  independent 
of  the  instantaneous  location  of  the  separation  point.  For  convenience,  such  a  contour  was 
chosen  to  coincide  with  the  normal  grid  line  such  that  the  roll  angle  of  its  intersection  with 
the  cylinder  surface  was  90°.  This  line  is  the  baseline  of  the  centered  velocity  profiles  in 
both  graphs  of  Figure  2. 

The  tangential  vorticity  flux  across  this  line  was  numerically  integrated  at  each  time 
step  between  r  =  D/2  and  the  outward  boundary  of  the  computational  domain.  By 
convention,  it  will  be  denoted  JOu0dr  and  refened  to  as  the  "vorticity  flux"  throughout  the 
remainder  of  this  manuscript. 

Typical  indicial  responses  are  shown  in  Figure  5  for  the  separation  point  and  the 
vorticity  flux,  in  the  case  of  a  1%  velocity  increase  at  Mc  =  0.25.  In  both  graphs,  t  =  0  is  the 
time  at  which  the  cylinder  was  set  in  motion.  The  two  graphs  are  plotted  on  the  same  time 
scale,  so  that  the  time  responses  A0$  and  AJOu0dr  can  be  easily  compared.  It  is  evident  that 
the  separation  roll  angle  reacts  to  the  instantaneous  change  in  velocity  with  a  relatively 
short  time  constant:  after  undergoing  a  sharp  positive  spike  (instantaneous  reattachment), 
the  separation  point  relaxes  back  to  a  position  slightly  upstream  of  the  initial  separation 
angle.  The  time  constant  associated  with  the  relaxation  process  is  seen  to  be  quite  small, 
since  0S  approaches  its  asymptotic  value  within  a  fraction  of  a  unit  convection  time. 

The  vorticity  flux,  on  the  other  hand,  experiences  a  significantly  longer  transient 
behavior,  and  exhibits  a  broad  overshoot  after  the  initial  step.  As  can  be  anticipated  from 
Equation  (4),  0$  and  JOu0dr  will  therefore  display  significantly  different  responses  to 
arbitrary  stimuli. 

In  order  to  establish  the  validity  bounds  of  the  indicial  theory  described  in  Section  2, 
the  test  disturbance  £  =  AMC  is  taken  to  be  of  the  form  elk  fc,  where  i2  =  -1  and  k'  is  the 
reduced  angular  frequency  based  on  diameter  and  crossflow  velocity.  The  advantages  of 
using  a  sinusoidal  disturbance  include:  i)  ease  to  quantify  the  accuracy  of  the  prediction 
through  a  decomposition  into  phase  and  amplitude  responses,  and  ii)  establishing  frequency 
bounds  for  the  application  range  of  Equation  (4).  One  of  the  objectives  of  the  present  study 
is  to  determine  whether  indicial  theory  can  be  used  to  predict  the  unsteady  motion  of  the 
separation  angle  and  the  unsteady  vorticity  flux.  Such  predictions  ought  to  be  possible  if 
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Equations  (12)  and  (13)  are  the  indicial  method  prediction  of  phase  delays  and  attenuation 
factors,  as  a  function  of  the  angular  reduced  frequency  k'. 


The  accuracy  of  this  prediction  was  tested  against  the  solution  of  actual  integrations  of 
the  Navier-Stokes  equations  at  specific  frequencies.  These  frequencies  were  typically  chosen 
to  be  at  critical  points  (i.e.  extrema,  inflection  points,  etc...)  of  the  analytically  predicted 
response  curves,  or  based  on  the  range  of  reduced  frequencies  of  practical  interest.  The 
most  severe  unsteady  effects  (highest  frequencies)  occur  in  situations  of  missile  flutter, 
rather  than  maneuver.  In  this  worst  case  scenario  the  range  of  amplified  frequencies9  is 
5Hz  <  f  <  30Hz.  Let  k  =  2/rfL/U^  be  the  reduced  angular  frequency  based  on  body  length 
and  freestream  velocity,  then  at  an  angle  of  attack  a,  the  reduced  frequency  seen  in  the 
crossflow  plane  (i.e.  based  on  crossflow  velocity  and  diameter)  is: 
k'  =  27TfD/Uc  =  27TDa^  (sincr)"^  f/M^,.  The  reduced  frequency  k'  is  therefore  maximized 
with  maximum  frequency,  f,  maximum  diameter  D,  minimum  freestream  Mach  number, 
and  minimum  angle  of  attack.  Since  for  typical  missile  aspect  ratios,  steady  separation  does 
not  occur  until  a  >  10°,  one  may  deduce  that  a  value  k'  -  10  represents  the  "flutter 
boundary".  For  this  reason,  direct  numerical  verification  of  frequency  effects  were 
constrained  to  the  range  k'  <  20. 

A  comparison  of  analytically  predicted  and  numerically  computed  phase  delays  and 
attenuation  factors  is  presented  in  Figure  6,  for  the  vorticity  flux.  The  solid  line  represents 
the  prediction  based  on  the  convolution  (Equation  (4))  with  the  indicial  response  of  the 
vorticity  flux  (Figure  5),  corresponding  to  a  1%  step  in  Mach  number.  The  equivalent 
prediction  corresponding  to  a  10%  change  in  velocity  (computed  to  evaluate  linearity 
bounds)  is  also  displayed  for  comparison  (dashed  line).  In  both  cases,  the  predictions  are 
based  on  an  "averaged"  indicial  response  (i.e.  the  average  of  the  response  to  a  positive  step 
on  the  one  hand  and  the  negative  of  the  response  to  a  negative  step  on  the  other).  In  the 
present  case,  the  "positive"  and  "negative"  responses  were  fairly  symmetric,  so  that  the 
averaging  process  is  not  fundamental  to  the  results,  but  rather  was  performed  for 
consistency  of  the  data  processing  with  cases  of  more  pronounced  asymmetry. 

As  may  be  seen  from  Figure  6  the  results  of  the  numerical  integrations  (symbols) 
follow  closely  the  analytical  curves,  with  an  overall  better  agreement  with  the  10% 
amplitude  case.  This  may  be  due  to  the  higher  resolution  of  the  large  amplitude  case.  In 
any  event,  a  significant  (50%)  overshoot  of  the  vorticity  flux  amplitude  is  predicted  and 
observed  at  a  frequency  k'  =  2,  which  is  well  within  the  range  of  observed  flutter 
frequencies.  A  phase  lag  of  approximately  30°  is  reached  around  k'  =  5.  Both  of  these 
dynamic  effects  may  have  profound  consequences  for  flow  prediction  using  quasi-steady 
implementations  of  separation  criteria. 

Unlike  the  vorticity  flux,  the  indicial  responses  (positive  and  negative)  of  the 
separation  angle  0S  become  asymmetric  at  large  amplitudes,  as  illustrated  in  Figure  7.  As 
will  be  shown  shortly,  this  asymmetry  (which  reflects  nonlinearity)  is  associated  with  a 
severe  deterioration  of  the  prediction  accuracy  of  the  indicial  method.  The  positive  step 
case  (cylinder  moving  against  the  flow  at  10%  of  the  freestream  velocity)  is  characterized  by 
the  instantaneous  creation  of  vorticity  at  the  body  surface,  which  translates  intc  a  negative 
vorticity  gradient  at  the  wall. 
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Hence  the  flow  reattaches  instantly  (according  to  the  detection  criterion  presented 
above)  at  time  t  =  0.  This  spike  in  A 9S  is  followed  by  a  long  transient  (of  the  order  of  one 
unit  convective  time),  corresponding  to  a  period  of  upstream-moving  separation.  When  the 
cylinder  is  impulsively  translated  with  the  flow  (dashed  line),  the  boundary  layer  instantly 
detaches  (negative  spike),  and  slowly  recovers  to  a  separation  angle  located  downstream  of 
the  original  location,  a  comparison  with  the  small  amplitude  step  response  of  Figure  5 
indicates  that  the  transient  time  increases  with  amplitude.  At  small  amplitudes  (1%  or  less), 
the  indicial  responses  in  0S  were  observed  to  be  symmetric  within  a  few  percent.  At  larger 
amplitudes,  the  asymmetry  seems  to  stem  from  the  fact  that  the  speed  at  which  the 
separation  point  travels  along  the  surface  depends  on  whether  the  separation  is  "upstream- 
moving"  or  "downstream-moving". 

The  predictive  accuracy  of  the  indicial  method  for  the  separation  point  location  was 
evaluated  against  selected  numerical  experiments,  similarly  to  the  vorticity  flux.  The  results 
of  these  comparisons  is  shown  in  Figure  8.  As  expected  from  the  above  discussion,  the 
computations  at  the  low  (one  percent)  amplitude  are  in  good  agreement  with  the  solid- 
curve  prediction.  At  the  low  amplitude  the  shape  of  the  indicial  response  was  observed  to 
be  independent  of  the  direction  of  the  step,  which  is  consistent  with  the  idea  of  local 
linearity  which  underlies  the  applicability  of  indicial  theory  as  outline  in  Section  2. 

For  the  large  amplitude  oscillations  (i.e.  peak  cylinder  velocity  equal  to  10%  of 
freestream),  a  significant  amount  of  harmonic  distortion  of  the  separation  angle  was 
observed,  as  shown  in  the  inset  of  Figure  8.  The  harmonic  distortion  is  indicative  of  the 
nonlinearities  which  are  also  reflected  in  the  asymmetry  illustrated  in  Figure  7.  Because  the 
numerically  computed  time-series  0s(t)  covered  too  few  cycles  of  the  fundamental  (forcing) 
oscillation,  adequate  narrow-band  frequency  filtering  was  not  attempted.  Instead,  the 
amplitude  of  the  oscillation  was  simply  defined  as  the  half  of  the  peak  to  peak  amplitude, 
while  the  phase  was  calculated  from  the  average  of  temporal  offsets  of  minima  and  maxima 
between  forcing  and  response  signals.  Due  to  the  severe  asymmetry  of  the  waveforms,  the 
uncertainty  on  the  phase  delays  thus  measured  was  considerable,  as  indicated  by  the 
uncertainty  bars  in  Figure  8.  Despite  the  uncertainty,  it  may  be  seen  that  the  estimated 
phase  delay  remains  close  to  the  predicted  curve  for  k'  <  10.  At  the  highest  frequency,  the 
numerical  solution  appears  to  follow  more  closely  the  small  amplitude  prediction.  For  the 
attenuation  factor,  A(k')/ A(0),  the  prediction  remained  within  25%  of  measurement  at  all 
frequencies. 

In  spite  of  these  inaccuracies,  several  important  features  may  be  extracted  from  the 
results  of  Figure  8.  First,  there  is  a  very  rapid  phase  adjustment  from  0°  to  -90°  occurring 
before  k'  =  1.  As  will  be  demonstrated  in  Appendix  B,  this  adjustment  can  be  theoretically 
predicted  based  on  a  simple  model  for. the  indicial  response  involving  the  superposition  of  a 
Heaviside  step  function  with  an  exponential  relaxation  curve.  It  may  be  shown  in 
particular,  that  the  width  of  this  frequency  window  is  inversely  proportional  to  the  time 
constant  of  the  relaxation  phase  and  to  the  square  root  of  the  amplitude  of  the  initial  pulse. 
The  narrowing  of  the  adjust menl  window  with  increasing  amplitude  is  visible  in  Figure  8 
Secondly,  the  attenuation  factor  is  seen  to  be  larger  than  unity  at  all  frequencies,  indicating 
an  amplification  of  fluctuations  of  the  separation  angle.  As  shown  in  Figure  8,  a  tenfold 
amplification  factor  can  be  obtained  at  relatively  low  reduced  frequencies,  underlining  once 
again  the  importance  of  unsteady  effects  on  the  dynamics  of  flow  separation. 
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As  a  final  test  of  the  applicability  of  indicial  theory  to  boundary  layer  separation,  an 
investigation  was  also  made  into  the  prediction  of  oscillatory  drag  for  viscous  flow  around 
a  two-dimensional  cylinder.  The  goal  of  this  investigation  was  to  establish  whether 
concepts  of  indicial  theory  would  also  be  applied  to  predict  oscillatory  lift  on  a  missile-like 
body  undergoing  flutter. 

Typical  indicial  responses  are  shown  in  Figure  9  for  viscous  and  pressure  drag 
components,  as  well  as  for  the  total  drag.  It  may  be  seen  that  although  it  is  the  viscous  drag 
which  is  responsible  for  the  ultimate  decrease  in  total  drag  coefficient  with  increasing 
velocity,  the  overall  character  of  the  transient  response  is  imposed  by  the  behavior  of  the 
pressure  drag.  While  the  viscous  drag  indicial  response  is  characterized  by  an  extremely 
narrow  initial  pulse,  both  pressure  and  total  drags  qualitatively  exhibit  the  same  features  as 
the  indicial  response  of  the  separation  angle.  The  existence  of  an  initial  adjustment  window 
for  the  phase  response  should  therefore  be  expected,  as  well  as  the  rapid  amplification  of 
fluctuations  with  reduced  frequency  (a  feature  also  predicted  by  the  model  discussed  in 
Appendix  B). 

This  is  numerically  verified  in  Figure  10.  In  the  top  graph  the  quasi-steady  phase  of  the 
pressure  drag  is  -180°,  reflecting  a  sign  inversion  that  was  applied  so  that  the  final  (quasi¬ 
steady)  value  of  all  drag  component  indicial  response  have  the  same  sign,  in  order  to 
facilitate  the  comparison.  At  all  tested  amplitudes  and  frequencies,  the  agreement  is  seen  to 
be  excellent. 

5.  Concluding  Remarks 

It  has  been  shown  for  viscous  flow  about  a  two-dimensional  cylinder,  that  several  key 
aspects  characterizing  the  time-dependent  behavior  of  boundary  layer  separation  can  be 
predicted  within  a  reasonable  degree  of  accuracy  over  a  large  range  of  frequencies.  The 
analytical  prediction  involved  only  a  convolution  integral  based  on  the  knowledge  of  the 
step  response  of  the  flowfield  to  a  small  perturbation.  The  perturbations  under 
consideration  were  changes  in  Mach  number  and  Reynolds  number.  The  test  functions  that 
were  used  to  evaluate  the  range  of  applicability  of  indicial  theory  to  the  separation  process 
were:  the  separation  angle,  the  vorticity  flux  across  the  separating  boundary  layer,  and  the 
drag  coefficient. 

Although  the  overall  agreement  between  the  indicial  method  prediction  and  the  results 
of  direct  numerical  simulations  was  clearly  superior  in  the  case  of  integral  quantities  such  as 
drag  and  vorticity  flux,  the  location  of  the  separation  point  was  found  to  be  satisfactorily 
predicted  within  the  range  of  frequencies  corresponding  to  missile  flutter.  The  breakdown 
of  the  method  appears  to  coincide  with  the  advent  of  nonlinearity,  which  is  first  conveyed 
by  an  asymmetric  response  to  positive  and  negative  step  inputs. 

Although  the  present  results  are  limited  to  low  Reynolds  number,  laminar  flow,  they 
establish  for  the  first  time  that  the  unsteady  characteristics  of  two-dimensional  flow 
separation  may  be  predicted  using  indicial  theory.  Whether  such  ideas  are  applicable  to 
high  Reynolds  number  turbulent  flow  remains  to  be  established. 
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APPENDIX  A 


PERTURBATIONS  IN  REYNOLDS  NUMBER 

It  was  shown  in  Section  3  that  the  process  of  impulsively  translating  the  entire 
computational  grid  is  that  closest  to  a  crossflow  plane  simulation  of  a  change  in  angle  of 
attack.  In  this  section,  the  concept  of  perturbing  a  separated  flow  solution  of  the  Navier- 
Stokes  equations  is  pursued  further  by  considering  the  effect  of  Reynolds  number  change. 
As  explained  earlier,  the  Reynolds  number  in  the  ARC2D  code  is  that  based  on  cylinder 
diameter  and  speed  of  sound.  Hence,  any  change  in  Re  =■  MOTRea  at  a  fixed  Mach  number 
should  be  interpreted  as  either  a  change  in  cylinder  diameter  or  (more  appropriately  in  the 
present  case)  as  a  change  in  kinematic  viscosity.  The  present  section  addresses  the  problem 
of  evaluating  the  limits  of  applicability  of  indicial  theory  when  a  solution  is  perturbed  in 
Reynolds  number. 

Figure  11  documents  the  change  in  separation  angle  0S  as  a  function  of  Reynolds 
number,  for  steady-solutions  of  the  Navier-Stokes  equations.  Thus,  the  Reynolds  number 
can  be  chosen  to  control  the  initial  angle  at  which  the  boundary  layer  separates. 

A  methodology  similar  to  that  described  in  Section  4  was  employed  here.  The  impulse 
response  of  the  separation  angle  0$  is  first  recorded.  The  responses  associated  with  positive 
and  negative  steps  in  Re  are  suitably  averaged  and  used  with  a  hypothetical  Reynolds 
number  oscillation  in  the  convolution  integral,  Equation  (4).  From  Equation  (11),  phase  and 
attenuation  frequency  responses  are  obtained  (Equations  (12)  and  (13)).  These  prediction 
curves  are  then  compared  to  the  results  of  direct  numerical  simulations  where  Re  is 
oscillated  sinusoidally  about  a  mean. 

Figure  12  presents  such  results  for  the  phase  delay  prediction  at  various  Reynolds 
numbers.  Since  Reynolds  number  and  separation  angle  are  quasi-statically  anticorrelated 
(see  Figure  11),  the  phase  difference  A 0  actually  involves  the  negative  of  the  Reynolds 
number  variation,  so  a  to  have  A0  =  0°  at  k'  =  0.  The  phase  prediction  is  seen  to  be  highly 
accurate  at  Re  =  40  and  Re  =  60.  In  these  two  graphs  the  dashed  line  indicates  a  simplified 
prediction  replacing  Equation  (11).  This  prediction  substitutes  an  exponential  least  square 
fit  to  the  actual  indicial  response,  i.e.  it  is  assumed  that  the  indicial  response  is  of  the  type 
0$(t)  =  0$(O)  +  (0$(«>)-0s(O))(l  -  e'bt).  It  is  then  easily  shown3  that  in  response  to  a  purely 
periodic  perturbation  e(t)  of  angular  reduced  frequency  k',  the  asymptotic  response  0s(t),  as 
t  ->  co,  will  exhibit  a  phase  lag  equation 

A0  =  -atan(k'/b)  (14) 


and  an  amplitude  attenuation  factor 


A(k')/A(0) 


(15) 


It  is  clear  that  the  simplified  prediction  of  Equation  (14)  becomes  at  the  higher  frequencies, 
therefore  stressing  the  importance  of  the  details  of  the  indicial  response. 

At  Re  =  90  (bottom  graph.  Figure  12),  the  prediction  is  quite  poor.  A  possible 
explanation  resides  in  the  reduced  resolution  in  the  determination  of  0$.  As  may  be  seen 
from  Figure  11,  the  decreased  sensitivity  d0s/dRe  at  Re  =  90  implies  total  angular  variations 
A0S  of  only  a  few  degrees,  which  is  of  the  order  of  the  grid  resolution  itself.  The  lack  of 
spatial  resolution  is  illustrated  in  Figure  13.  The  dashed  line  represents  the  average  of  the 
two  roll  angles  delimiting  the  computational  cell  in  which  the  normal  vorticity  gradient 
80/  0n  crosses  zero  at  any  point  in  time.  The  solid  lines  represent  the  actual  (positive  and 
negative)  indicial  resi.  ^nses,  obtained  by  linear  interpolation  of  00/  0n.  Finally,  the  dotted 
lines  represent  crude  exponential  fits  to  the  indicial  responses.  Figure  13  shows  that  a 
significant  amount  of  asymmetry  is  already  present  at  Re  =  40  but  does  not  seem  to 
significantly  impair  the  accuracy  of  the  prediction  (top  graph.  Figure  12).  At  Re  =  90,  a 
deterioration  of  both  the  angular  resolution  and  the  degree  of  symmetry  of  the  step 
responses  was  observed.  It  is  believed  that  these  factors  account  for  the  poor  agreement 
observed  in  Figure  12  (bottom  graph).  The  influence  of  the  asymmetry  between  positive 
and  negative  indicial  responses  is  also  shown  in  this  Figure. 

Figure  14  depicts  similar  comparisons  for  the  attenuation  factor  of  the  oscillatory 
separation  angle  at  Re  =  40,  Re  =  60,. and  Re  =  90.  The  validity  of  the  simplified  predictions 
based  on  an  exponential  fit  of  the  indicial  response  applies  to  a  smaller  frequency  range  than 
for  the  case  of  phase  lag  prediction.  The  numerically  computed  results  are  in  qualitative 
agreement  with  the  full  analytical  prediction.  The  apparently  large  errors  (-20%)  may  be 
attributable  to  the  lack  of  spatial  resolution  for  all  cases  involving  Reynolds  number 
perturbations. 

Figures  15  and  16  depict  the  phase  and  attenuation  results  at  a  fixed  initial  Reynolds 
number,  for  crossflow  Mach  numbers  of  Mc  =  0.1,  Mc  =  0.25,  and  Mc  =  0.4.  At  the  lowest 
Mach  numbers,  the  phase  lag  prediction  are  in  excellent  agreement  with  the  results  of  the 
direct  numerical  simulations.  However,  in  the  case  Mc  =  0.4,  both  predicted  phase  lags  and 
attenuation  factors  strongly  deviate  from  the  numerical  results.  As  in  the  higher  Reynolds 
number  case,  the  importance  of  the  asymmetry  between  positive  and  negative  step  response 
is  indicated  by  the  corresponding  individual  predictions  in  Figures  15  and  16  (bottom).  In 
this  case,  the  origin  of  the  asymmetrical  responses  is  not  known,  although  is  may  be 
conjectured  that  compressibility  effects  may  have  come  into  play,  since  shocks  are  known  to 
appear  at  freestream  Mach  numbers  as  low  as  =0.55.  In  any  event,  the  lack  of  spatial 
resolution  may  also  have  played  an  important  role,  since  the  total  variation  in  static 
separation  angle  0S  between  Mc  =  0.1  and  Mc  =  0.4  at  Re  =  60  was  less  than  0.7°,  in 
comparison  with  a  total  spread  of  approximately  9°  between  Re  =  40  and  Re  =  90  at 
Mc  =  0.1. 

Finally  it  is  worthwhile  noting  that  the  phase  response  of  the  separation  angle  shows 
only  a  weak  dependence  on  Reynolds  number,  but  a  relatively  strong  dependence  on  the 
Mach  number.  This  can  be  seen  by  comparing  the  evolution  of  phase  lag  distributions 
between  Figures  12  and  15.  Furthermore,  the  time  constants  of  the  exponential  fits  were 
found  to  be  numerically  close  in  the  case  of  the  Reynolds  number  variation,  but  exhibited  a 
definitive  trend  with  Mach  number,  as  shown  in  Figure  17. 
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To  test  the  ability  of  the  indicial  method  to  predict  the  dynamic  characteristics  of  the 
vorticity  flux  with  oscillatory  Reynolds  number,  a  sample  calculation  was  performed  in  the 
vicinity  of  the  frequency  which  indicated  theoretically  (i.e.  according  to  Equation  (11))  a 
150%  amplification  in  fluctuating  amplitude.  As  shown  in  Figure  18,  the  predicted  behavior 
faithfully  reproduced  the  results  of  the  numerical  simulation.  The  presence  of  an  overshoot 
at  the  same  frequency  of  k'  =  2.5  was  previously  noticed  in  Section  4  under  the  same 
nominal  conditions  of  Mach  and  Reynolds  number,  but  for  totally  different  types  of 
perturbations. 
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APPENDIX  B 


ANALYSIS  OF  A  MODEL  INDICIAL  RESPONSE 


Let  £(t)  designate  the  indicial  response  of  a  function  //(t)  in  response  to  a  perturbation 
£,  then,  under  the  conditions  outlined  in  Section  2,  Equation  (4)  must  hold: 


t 


/i  (t)  =  f(t)e(0)  - 

o 


{(T)—  £  (t  -  T)dT 


(16) 


Let  H  designate  the  Heaviside  step  function.  If  the  following  form  is  assumed  for  the 
indicial  response: 


£<t)  =  [(l  - 


e  bT)  +  fjH(T) 


(17) 


then  Equation  (17)  is  seen  to  represent  a  step  function  of  magnitude  to  which  is 
superimposed  an  exponential  relaxation  function  of  unit  amplitude,  the  ratio  (£o/0+£o)) 
therefore,  measures  the  ratio  of  the  initial  value  (or  "spike",  if  <  0)  to  the  final  value  (i.e. 
after  all  transients  have  decayed). 


Substituting  Equation  (17)  into  Equation  (16)  and  assuming  the  perturbation  e  to  be  of 
purely  oscillatory  form  (e(t)  =  eiwt,  i2  =  -1)  yields  for  large  times: 


iwt  r„  . 

n(t)  =  e  [l  +  £  " 


Hence  the  phase  response  $  is  given  by: 

( p  =  -atan 

where 


iw(b  -  iw) 

2  TT 

w  +  b 


•] 


w 


b  +  c 


c  = 


2  ,2 

w  +  b 


'  b 


while  the  attenuation  factor  can  be  expressed  as 


(18) 


(19) 


(20) 


A(w)  |n  . 

1 

A(0)  'Ii  " 

1  + 

w 

b 

2-  1 


(21) 


where 
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2 


1  + 


(22) 


At  finite  frequency.  Equation  (21)  becomes  singular  at  =  -1.  This  situation  arises  when 
the  final  value  of  fi  equals  its  amplitude  immediately  before  the  initiation  of  the  step.  This 
does  not  typically  occur  in  practice.  However,  it  is  useful  to  examine  the  behavior  of  the 
amplitude  (attenuation)  response  in  the  case  where  the  indicial  response  is  characterized  by 
a  large  initial  spike  followed  by  a  relaxation  phase  back  to  a  value  which  is  very  close  to  the 
initial  value.  This  type  of  indicial  response  was  observed  for  the  separation  angle  0S  and  the 
drag  coefficient  (Figures  7  and  9),  and  may  be  examined  in  detail  by  considering  first  order 
expansions  of  Equations  (19)  and  (21)  about  £0  =  -1. 

Let  £0  =  -1  +  5,  where  I  6  I  «1.  It  may  be  shown  that 


A(w) 

A  (0)  S->  0 


(23) 


Hence,  the  shape  of  the  amplitude  response  curve  is  given  by  the  square  root  in  (23),  while 
at  any  given  frequency,  the  response  A(w)/A(0)  gets  amplified  by  a  factor  6"1,  i.e.  the 
magnitude  of  the  spike  of  the  indicial  response,  relative  to  the  final  value.  This  behavior  is 
exemplified  by  the  amplitude  response  curves  associated  with  pressure  and  total  drag 
(Figure  9). 


A  similar  analysis  of  the  phase  response  (Equation  (19))  reveals  that  the  phase 
asymptotes  to  -180°  for  large  frequencies  (a  fact  corroborated  by  the  present  results.  Figure 
10,  top  graph),  and  that  the  initial  rate  of  change  of  the  phase  lag  <j>  with  respect  to  w  is 


d£  _ > 

dw  S->0 

w=0 


1_ 
b  8 


(24) 


It  can  be  shown  that  for  <5  fixed  and  w  0,  then  <p  -» 0.  However,  in  the  limit  where 
0  <  w  «  1  is  fixed  and  6  *■»  0,  the  phase  asymptotes  instead  to  a  value  of  -7t/2.  Since  6  is 
never  zero  in  practice,  it  is  the  zero  limit  which  is  observed.  Nevertheless,  $ (w)  may  be 
viewed  as  the  composite  of  two  solutions  which  match  at  a  frequency  w*  such  that  0(w*)  =  - 
n/2.  Therefore  w  represents  the  frequency  window  of  rapid  adjustment  of  the  phase 
between  $  =  0  and  $  =  -n/2.  The  frequency  w*  at  which  tan$  becomes  singular  is  given  by 
the  estimate 


w*  TS*  hU  (25) 

Hence  for  a  fixed  relaxation  constant  b"1,  the  width  w*  of  the  region  of  rapid  phase 
adjustment  decreases  like  the  square  root  of  the  spike  amplitude. 


-16- 


A  wealth  of  additional  information  can  be  extracted  from  such  simple  models  as  that 
expressed  in  Equation  (17).  In  particular,  the  analysis  can  be  easily  extended  to  the  case 
£  >  0,  corresponding  to  the  initial  transient  of  the  vorticity  flux  indicial  response  (Figure  5). 

The  analysis  of  the  model  indicial  response  (17)  with  £0  >  0  did  not  yield  the  observed 
overshoot  in  the  attenuation  curve  of  the  vorticity  flux.  Since  such  a  model  is  unable  to 
represent  the  local  maximum  exhibited  by  the  indicial  response  itself,  it  is  likely  that  this 
maximum  is  the  cause  of  the  overshoot.  In  any  event,  the  sensitivity  of  the  results  to  the 
magnitude  £0  was  found  to  be  small  for  £0  >  0,  but  very  significant  for  £0  <  0  (i.e.  in  the 
presence  of  an  initial  spike). 

Since  the  magnitude  of  the  initial  spike  appears  to  condition  strongly  the  phase  and 
amplitude  responses  for  9S  and  CD,  an  attempt  at  determining  its  scaling  properties  with  • 
respect  to  the  step  amplitude  was  performed.  Using  the  perturbation  technique  described 
in  Section  3,  the  Mach  number  was  stepped  up  from  an  initial  value  Mc  =  0.25  with  various 
amplitudes  covering  five  orders  of  magnitude,  from  AM/Mj  =  0.001%  to  AM/M.  =  100%. 
The  results  are  presented  on  a  logarithmic  scale  in  Figure  19.  In  both  graphs,  the  subscript  i 
designates  the  initial  value  (i.e.  immediately  before  the  cylinder  was  set  in  motion).  The 
round  symbols  represent  the  peak  value  of  the  indicial  responses  (at  t  =  0+)  while  the  square 
symbols  correspond  to  the  post-transient  value  (t  +°°).  The  solid  curves  are  lines  of  unit 
slope  which  were  added  as  an  aid  in  determining  the  range  of  linearity  of  the  data. 

Both  base  and  peak  values  for  the  indicial  response  of  the  total  drag  coefficient  behave 
very  linearly  with  forcing  amplitude.  Therefore,  the  ratio  (fo/0  +  £0))  is  a  constant. 
Consequently,  if  the  relaxation  constant,  b,  is  a  constant,  the  phase  and  amplitude 
characteristics  of  CD  should  be  independent  of  Mach  number  change  over  the  investigated 
range.  This  result  is  essentially  in  agreement  with  the  data  of  Figure  10.  In  the  case  of  the 
separation  angle  0$,  an  approximate  linearity  range  exists  below  AM/Mi  =  1%,  followed  by 
a  saturation  of  the  response.  Therefore,  it  is  expected  that  the  phase  and  attenuation 
characteristics  should  exhibit  a  larger  dependence  on  the  initial  perturbation  magnitude,  a 
•fact  also  corroborated  by  the  numerical  experiments  reported  in  this  manuscript. 
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Fig.  1.  Close-Up  View  of  Two-Dimensional  Grid  Used  for  Navier  Stokes  Calculations 
of  Unsteady  Flow  Around  a  Circular  Cylinder. 
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Fig.  5.  Typical  Indicial  Responses  for  the  Separation  Point  Qs  (Top)  and  the  Vorticity 
Flux  (Bottom),  for  Flow  Around  a  Circular  Cylinder  Subjected  to  a  1  %  Change 
in  Velocity  (M  =  0.25,  Re  =  60). 
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•Pig.  6.  Comparisons  Between  Theoretical  Predictions  (Based  on  the  Indicia! 
Responses)  and  hTumerical  Experiments  for  the  Vorticity  Flux.  Top  Graph: 
Phase  Delays  Bottom  Graph:  Amplitude  Attenuation  Factors.  (M  &  0.25 
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Fig.  8.  Comparisons  Between  Theoretical  Predictions  (Based  on  the  Indicial 
Responses)  and  Numerical  Experiments  for  the  Separation  Angle.  Top  Graph: 
Phase  Delays.  Bottom  Graph:  Amplitude  Attenuation  Factors.  (Mc  =  0.25, 
Re  =  60).  The  Inset  Illustrates  the  Harmonic  Distortion  which  Occurs  when 
Positive  and  Negative  Indicial  Responses  are  Not  Symmetric  (From  Direct 
Navier  Stokes  Simulation). 


Indicial  Responses  of  Viscous,  Pressure  and  Total  Drag  for  a  1%  Change  in 
Velocity  (Mc  =  0.25,  Re  =  60). 
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Fig.  10.  Comparisons  Between  Theoretical  Predictions  (Based  on  the  Indicial 
Responses)  and  Numerical  Experiments  for  Viscous,  Pressure,  and  Total  Drag. 
Top  Graph:  Phase  Delays.  Bottom  Graph:  Amplitude  Attenuation  Factors. 
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Fig.  12.  Evolution  of  Prediction  Accuracy  of  Phase  Delays  of  0S  with  Reynolds 
Number,  fqr  a  10%  Change  in  Reynolds  Number  at  Mc  =  0.1. 
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Evolution  of  Prediction  Accuracy  of  Attenuation  Factors  of  0S  with  Reynolds 
Number,  for  a  10%  Change  in  Reynolds  Number  at  Mc  =  0.1. 


PHASE 


I - i _ 1 - 1 _ ! _ i _ 1 _ i _ &  _ I _  ■  _ F 

0  5  10  15  20  25  30  35 

k'  =2TTfD/U, 


30 


o> 

<D 

CS' 

-30 

$ 

-60 

-90 

-120 

-  *predictlon  bhsed  on  convolution  wit#»  avenged  lAdlcial  response 

•  •  •  •  prediction  based  on  convolution  with  exponential  fit  to  the  indicial  resp- 
— •  prediction  based  on  convolution  with  negative-step  indicial  response 
—  prediction  based  on  convolution  with  positive-step  indicial  response  _ 
numerical  experiment 
-  ' 

1 

\ 

\ 

\ 

\ 

. . . .  „ 

1 

l 

V\\ 

— 1 — 1 — 

-  II 
o 
■  k 

I 

i 

i 

,  Ia 

\ 

\  f\  - 
!  ,  \;  \ 

O  5  10  15 

20 

25  30  3 

k'  =  2tt  fD/U 


Fig.  15.  Evolution  of  Prediction  Accuracy  of  Phase  Delays  of  0S  with  Cross-Flow 
Mach  Number,  for.  a  10%  Change  in  Reynolds  Number  at  Re  =  60. 


A(k')/A(0)  A(k')/A(0)  A(k-)/A(0) 


-  prediction  based  on  convolution  with  averaged  indicial  response 
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Fig.  16.  Evolution  of  Prediction  Accuracy  of  Attenuation  Factors  of  Q$  with  Cross- 
Flow  Mach  Number,  for  a  10%  Change  in  Reynolds  Number  at  Re  =  60. 
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Fig.  19.  Determination  of  the  Linear  Range,  with  respect  to  Mach  Number  Change,  of 
Peak  and  Post-Transient  Values  of  the  Indicial  Responses  for:  Total  Drag  (Top) 
and  Separation  Angle  (Bottom),  (Mc  =  0.25,  Re  =  60). 
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Abstract 


1.  Introduction 


The  Transonic  Small  Disturbance  equation  was 
supplemented  with  a  transport  equation  for  the 
streamwise  vorticity  and  a  vector  potential  equation  to 
predict  vortex  effects  over  missile  configurations.  The 
flow  separation  phenomenon  was  modeled  using  normal 
vorticity  jets  placed  along  the  separation  line.  The 
strength  and  location  of  the  separating  vorticity  was 
determined  from  empirical  criteria.  Time-accurate 
calculations  performed  using  a  modified  version  of  the 
CAP-TSD  code  in  subsonic,  transonic,  and  supersonic 
flow  suggest  that  it  is  possible  to  compute  realistic  angle 
of  attack  configurations  using  CAP-TSD,  thus  showing 
considerable  potential  for  aeroelastic  applications  and 
unsteady  aerodynamics. 

Nomenclature 


Symbols: 

A  Rotational  vector  potential 

CN  Normal  force  coefficient  (based  on  body  cross- 
sectional  area) 

Cp  Pressure  coefficient 

D  Body  diameter 

f  Frequency 

Im(}  Imaginary  part 

k  Angular  reduced  frequency  based  on  length  and 
freestream  velocity  CnfL/U'J 
k'  Angular  reduced  frequency  based  on  diameter 
and  crossflow  velocity  (PjrfD/liy 
L  Body  length 

M  Mach  number 

p  Pressure 

q  Total  velocity  (u^v^+w2)1  /2 

s  Curvilinear  coordinate  along  body  radius 

Re(}  Real  part 

t  Time 

u  Streamwise  velocity 

v  Spanwise  velocity 

w  Normal  velocity 

x  Streamwise  coordinate 

y  Spanwise  coordinate 

z  Normal  coordinate 

a  Angle  of  Attack 

T  Circulation 

0  Roll  angle 

(f>  Irrotational  potential 

p  Fluid  density 

O  Streamwise  vorticity 

Subscripts: 
c  Crossflow 

le  Leading  edge 

n  Normal 

te  Trailing  edge 

«  Freestream 
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In  recent  years  it  has  been  possible  to  predict  the 
unsteady  transonic  flow  around  a  wing,  especially  those 
typical  of  commercial  aircraft,  in  a  fairly  efficient 
manner.  Frequently,  the  computer  codes  that  are  used 
are  based  on  potential  theory  and  are  considerably  faster 
than  a  corresponding  calculation  using  the  Euler 
equations.  Because  these  methods  have  been  developed 
for  aircraft,  they  are  not  really  applicable  to  missile  flows 
where  the  effects  of  vorticity  due  to  flow  separation  are 
important;  the  potential  approximation  cannot  predict 
the  effects  of  vorticity  other  than  by  representing  a 
vortex  wake  by  an  infinitely  thin  sheet  which  is  excluded 
from  the  computational  domain.  This  model  is 
complicated  and  may  not  be  a  viable  option  for  routine 
calculations  around  real  aircraft  or  missile  configurations 
because  the  geometry  of  the  vortex  sheet  can  get  quite 
complex. 

For  steady  missile  flow  a  variety  of  prediction 
methods  are  available,1  ranging  from  panel  methods 
with  the  addition  of  nonlinear  vortex  dynamics  to  the 
Euler  equations  or  Navier-Stokes  equations.  It  would  be 
ideal  to  use  the  Navier-Stokes  equations  to  model 
unsteady  transonic  flow  around  missiles,  but  there  are 
several  difficulties  with  this  approach.  The  most  obvious 
difficulty  is  the  computer  time  required,  which  is  several 
orders  of  magnitude  greater  than  that  required  for  a 
potential  calculation.  Even  if  a  dedicated  supercomputer 
were  available  for  such  calculations,  the  computation  of 
the  necessary  flow  separation  might  be  inaccurate 
because  of  the  inherent  inaccuracy  in  many  present 
turbulence  models.  The  next  best  approach  would  be  to 
use  the  time  dependent  Euler  equations  which  require 
some  empiricism  to  initiate  separation  and  compute  the 
shed  vorticity.  If  the  separation  line  and  shed  vorticity 
could  be  predicted,  then  the  Euler  equations  would  be  a 
viable  model  since  they  will  model  the  transport  of  this 
vorticity  reasonably  accurately.  However,  since  missiles 
are  slender  some  further  approximations  can  be  made. 

This  paper  is  concerned  with  developing  a  method 
of  predicting  the  unsteady  pressure  distribution  on 
missile-like  bodies  at  transonic  and  supersonic  speeds. 
The  concepts  developed  in  this  work  are  extensions  of 
earlier  analyses2'3  for  steady  and  unsteady  flow.  The 
approach  is  to  make  as  much  use  of  the  existing 
technology  as  possible.  The  final  goal  is  a  computer  code 
capable  of  predicting  the  effects  of  three-dimensional 
unsteady  separation  in  transonic  flqw  about  missiles  for 
use  in  aeroelastic  calculations  or  during  maneuver.  One 
important  aspect  of  this  work  involves  the  application,  cf 
indicial  theory  to  predict  the  time-dependent  behavior  of 
boundary  layer  separation.  The  present  paper  is  an 
account  of  the  final  phase  of  the  work,  namely  the 
refinement  and  exploitation  of  the  unsteady  flow  variant 
of  this  method  for  prediction  of  three-dimensional 
unsteady  separated  flow. 
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2.  Basic  Equations 


3.  Computer  Code  and  Numerical  Methods 


The  present  approach  is  based  on  experience 
gained  in  steady  subsonic  and  supersonic  flow 
predictions  for  missiles,  in  particular,  the  fact  that  if  a 
separation  line  and  the  strength  of  the  vorticity 
introduced  at  that  line  can  be  estimated  empirically,  the 
governing  equations  (such  as  the  Euler  equations)  will 
represent  the  transport  of  this  vorticity  accurately. 
Because  the  computational  time  for  the  Euler  equations 
for  an  unsteady  calculation  is  considerably  greater  than 
that  for  a  potential  calculation,  a  simplified  model  is 
used. 


By  making  the  small  disturbance  approximation,  p  . 
can  be  expanded  in  terms  of  disturbance  velocities,  and 
substituted  into  the  conservation  of  mass  equation  to 
give  the  following  three-dimensional  unsteady  potential 
equation: 
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The  basic  equations  are  an  extension  of  those 
derived  by  Klopfer  and  Nixon4  for  steady  flow  and 
constitute  a  subset  of  the  Euler  equations.  For  slender 
bodies,  the  five  Euler  equations  are  reduced  to  an 
equation  similar  to  the  three-dimensional  unsteady 
potential  equation,  a  three-dimensional  vorticity 
transport  equation,  and  a  two-dimensional  Poisson 
equation.  In  addition  to  having  to  solve  a  reduced  set  of 
equations,  a  significant  advantage  is  that  one  of  the 
equations  is  almost  identical  to  the  potential  equation  for 
wnich  there  are  well  tested  computer  codes.  This  allows 
the  development  of  a  prediction  method  based 
considerably  on  proved  technology.  The  reader  is 
referred  to  Ref.  3  for  the  derivation  of  these  equations,  as 
well  as  a  description  of  their  numerical  implementation. 

The  basic  equations  are  summarized  as  follows: 
Vorticity  Transport 

flt  +  ^  +  Cvnjy +(wO)2==0  (1) 

Rotational  Crossflow 

Ayy  +  =  -ft  (2) 

Mass  Conservation 

pt  +  (pu)x  +  (pv)y  +  (pw)2  =  0  (3) 

where  p  is  the  fluid  density,  u,  v,  and  w  designate  the 
streamwise,  spanwise,  and  normal  velocity  components 
(i.e:,  in  the  coordinate  directions  x,  y  and  z,  respectively), 
and  A  is  the  vector  potential  for  the  rotational  velocity. 
The  Helmholz  decomposition  of  the  flowfield  into 
potential  and  rotational  components  is  given  by. 

u  ss  1  +  ;  v  sr  +  Az ;  w  =  $2  -  Ay  (4) 

where  is  the  irrotational  scalar  potential.  Finally,  the 
density  is  related  to  the  velocity  components  through  the 
isentropic  ideal  gas  relation: 
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where  Mtt  designates  the  freestream  Mach  number,  and 
Y  is  the  specific  heat  ratio. 

The  boundary  conditions  on  the  body  are  zero 
flow  through  the  body,  and  zero  vorticity  flux,  except  at 
the  specified  separation  line. 


f0=-M^t-2M^x 

f,  =  (1-Mi)^  -}<y+l)M2^  +^y-3)M2<#y+Az)2 
f2  =  #y+A2-(y-l)M^x«Sy  +A2) 

f3=^z'Ay 


(7) 


In  Ref.  3,  Nixon  et  al.  showed  that  because  of  the 
resemblance  between  Eq.  6  and  the  usual  Transonic 
Small  Disturbance  (TSD)  equation,  this  theory  could  be 
implemented  using  established  and  well-tested  potential 
flow  solvers,  such  as  the  Computational  Aeroelasticity 
Program  CAP-TSD,  aue  to  Batina  et  al.5  The  CAP-TSD 
code  is  capable  of. predicting  the  unsteady  transonic  and 
supersonic  flow  over  a  complete  aircraft  configuration 
including  fins,  stores,  and  pylons.  Consequently,  this 
code  was  used  as  a  basis  for  the  implementation  of  the 
above-mentioned  unsteady  flow  theory  (an 
implementation  of  the  steady  version  of  this  theory  may 
be  found  in  Ref.  4). 


The  present  paper  reports  sample  calculations  for  a 
model  three-dimensional  slender  body  with  horizontal 
and  vertical  fins.  The  results  illustrate  the  predictive 
capability  of  a  computational  tool  initially  based  on 
CAP-TSD,  but  which  has  extended  the  range  of 
applications  of  this  code  to  separated  vortical  flows, 
including  crucial  refinements  in  the  representation  of  the 
computational  fuselage. 

The  algorithm  used  in  CAP-TSD  is  approximate 
factorization  with  time  stepping  as  an  option;  this 
algorithm  is  retained  in  solving  Equation  (6)  and  the 
vorticity  transport  equation.  Equation  (1)*  Since 
Equation  (2)  does  not  have  explicit  time-dependence,  it 
can  be  solved  at  each  time  step  and  in  each  individual 
streamwise  plane  by  rising  the  successive  over-relaxation 
method.  The  time-marching  procedure  is  performed 
using  a  second  order  accurate  implicit  scheme.  For  the 
continuity  and  vorticity  transport  equations,  if  the 
Newton  iteration  option  is  used,  then  the  L^-norm  over 
space  of  the  incremental  change  in  the  solution  (A$  or 
AD)  approaches  zero  at  each  time  step,  and  the  solution  is 
time-accurate.  For  the  vector  potential  equation,  the 
solution  is  also  time-accurate,  since  at  each  iteration,  AA 
can  be  mad$  arbitrarily  close  to  zero.  A  more  complete 
description  of  the  algorithm  and  the  difference  equations 
which  have  been  used  to  modify  the  CAP-TSD  code  can 
be  found  in  Ref.  3., 
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4.  Boundary  Conditions 

a)  Body.  At  the  body#  the  /low  is  required  to  be 
tangential  to  the  surface,  and  a  normal  vorticity  jet  is 
prescribed  at  the  point  of  flow  separation. 

In  CAP-TSD,  the  fuselage  cross-section  is  assumed 
to  be  elliptic,  and  the  boundary  conditions  on  the  body 
surface  are  transferred  to  a  bounding  cartesian  grid  (the 
computational  fuselage)  by  using  ideas  from  slender 
body  theory.  For  angles  of  attack  typical  of  missile 
flight,  however,  the  original  treatment  was  inadequate 
and  an  alternative  scheme  was  developed.3  In  this 
modified  scheme,  the  assumption  of  slender  body  flow  is 
still  used,  in  particular,  the  consequence  that  only  .the 
crossflow  is  important.  The  boundary  conditions  on  the 
cartesian  computational  boundary  are  found  from  the 
'thickness  effect/  which  is  relatively  accurate,  and  by 
using  an  analytic  solution  for  the  crossflow  to  specify 
normal  velocities  v  or  w  on  the  computational  boundary. 
For  a  body  at  angle  of  attack  a  these  conditions,  in  a 
reference  frame- aligned  with  the  freestream, are 
expressed  as: 


+  Az=v  ;  *Z-Ay  =  w-sin(a)  (8) 

Since  the  standard  form3  of  the  TSD  Equation  (i.e., 
in  the  absence  of  swept  shocks)  reduces  to  Laplace's 
equation  in  the  crossflow  plane,  it  is  justified  to  use 
analytic  solutions  v  and  w  obtained  from  conformal 
mapping  transformations.  Such  an  approach  has  been 
previously  shown3  to  be  valid  for  reduced  frequencies, 
2;rfL/U w  ,  up  to  order  unity,  within  the  context  of 
slender  body  approximation  and  negligible 
compressibility  effects  (of  order  a2)  in  the  crossflow. 

For  a  body  with  elliptical  cross-section  and  thin 
lifting  surfaces  aligned  with  the  principal  axes  of  the 
ellipse,  this  analytic  solution  is  obtained  by  using 
successive  conformal  mapping  transformations.  The 
global  transformation  (noted  c  =  *(Z),  with  Z  =  -z  +  iy) 
maps  the  symmetry  plane  (y  =  0),  the  elliptical 
boundary,  and  the  horizontal  and  vertical  fins  into  a 
subset  of  the  real  axis  in  the  tr-plane.  The  transformation 
is  the  composite  of  the  following  mappings: 


i)  ellipse  of  semi-axes  b  and  a  into  the  unit  .circle: 


Z  +  (Z2-a2+b2)1/2 
a+b 


(9) 


ii)  unit  circle  into  a  flat  plate: 

n  “  £  +  j 

iii)  flat  plate  with  vertical  fence  of  height  rj  into  a 
flat  plate: 


In  the  (T-plane,  the  flow  is  that  associated  with  a 
point  vortex  in  a  crossflow  (see  Nixon  et  al.3).  The 
freestream  velocity  magnitude  in  <r-space  is 
(a+b)fjcsin(a)/2.  The  point  vortex  used  in  each 
stream  wise  plane  models  the  rotational  component  of 
the  flowfiela  by  being  given  a  strength  ro  equal  to  the 


numerically  integrated  vorticity  in  that  plane.  The 
location  of  the  vortex  in  Z-space  is  chosen  to  be  the 
centroid  of  vorticity.* 


Let  F  be  the  complex  potential  associated  with  the 
analytic  solution  in  the  transformed  plane;  the  spanwise 
and  normal  velocities  v  and  w  can  then  be  recovered 
according  to: 


V  - 


(12) 


and 


If  A  is  a  constant  along  the  boundary,  then 
Equations  8  are  used  to  supply  CAP-TSD  with  the 
required  Neumann  boundary  conditions  on  The 
inclusion  of  the  sin(a)  term  in  Equation  8  reflects  the  fact 
that  the  equations  of  motion  are  solved  in  a  reference 
frame  attached  to  the  body,  i.e.,  at  an  angle  of  attack  a 
from  the  freestream. 


In  the  original  CAP-TSD  code,  the  computational 
fuselage  is  a  rectangular  cylinder  which  encompasses  the 
entire  body.  Such  a  configuration  was  used  in 
computations  by  Nixon  et  al./3  for  a  body  with  a  2:1 
elliptical  cross-section.  For  such  bodies,  the  separation 
line  lies  close  to  the  90*  roll  angle,  and  the  shear  layers 
are  known  to  separate  almost  vertically.  There  are, 
however,  several  problems  with  this  representation  of 
the  fuselage,  in  particular,  the  fact  that  the  exact  location 
of  the  vorticity  jk  cannot  be  determined  accurately,  since 
one  must  account  for  the  transport  of  vorticity  within  the 
space  bounded  by  the  true  and  computational  body 
surfaces.  Furthermore,  the  vorticity  distribution  at  a  far 
away  boundary  cannot  be  represented  by  a  single 
outward  vorticity  jet  at  some  location  on  that  boundary. 
This  problem  becomes  especially  acute  with  low  aspect 
ratio  ellipses,  and  in  particular  in  the  case  of  circular 
cross-sections. 

The  difficulty  was  circumvented  by  modifying  the 
CAP-TSD  code  so  as  to  represent  the  computational 
fuselage  by  a  serrated-edged  cylinder  which  closely 
approximates  the  actual  solid  boundary.  The  obvious 
advantages  of  this  modification  are:  i)  an  accurate 
calculation  of  the  vorticity  transport  phenomenon  close 
to  the  separation  point;  ii)  a  lesser  dependence  on  the 
model  used  for  tne  transfer  of  boundary  conditions 
(Equations  12  and  13);  and  iii)  convergence  of  the  true 
and  computational  surfaces  towards  one  another,  with 
increasing  grid  resolution  in  the  crossflow  plane.  Most 
importantly,  this  representation  of  the  computational 
fuselage  can  retain  the  cartesian  grid  of  the  CAP-TSD 
code.  The  resulting  grid  is  shown  in  Fig.  1,  where  the 
cross-section  of  the  actual  body  has  been  outlined.  The 


_ _ _ 

•In  the  two-dimensional  point  vortex  model,  it  is  necessary  to  consider  the 
existence  of  a  'feeding  sheet'  along  the  branch  cut  because  ar/  iM). 
This  is  because  the  variable  circulation  results  in  a  af /at  term  in  the 
BemouiUi  equation,  hence  producing  a  pressure  jump  across  the  branch 
cut.  The  resulting  force  must  be  balanced  by  the  Joukowski  force  on 
the  vortex  core.*’  This  'force-free'  (but  not  moment-free)  vortex 
condition  leads  to  an  O.D.E.  for  the  vortex  position  and  circulation20  as 
a  function  of  time,  for  given  initial  conditions.  Since  r  Is  known  such 
an  approach  would  essentially  amount  to  a  duplication  of  the  vorticity 
transport  equation.  Therefore,  the  simple  centroid  model  was  retained 
in  the  present  boundary  conditions. 
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arrow  in  the  close-up  figure  indicates  the  roll  angle 
convention. 

For  the  vorticity  boundary  condition  at  the  body,  a 
flow  separation  condition  is  simulated  by  the  injection  of 
vorticity  into  the  flowfield  at  the  point  of  separation.  For 
steady  flow,  a  modified  version  of  the  Stratford 
criterion6  is  used  to  determine  the  separation  line: 

dc'  1/2.  (14> 

Cp  |R®3  x  10  I  *-0.35  sin  (a) 

where  Cp*  is  the  modified  pressure  coefficient7  and  s  is 
the  virtual  length  of  the  turbulent  boundary  layer,  as 
seen  in  the  crossflow  plane. 

Although  two-dimensional  Navier-Stokes 
calculations  of  unsteady  flow  around  a  circular  cylinder 
have  been  performed8  as  a  first  step  towards  the  goal  of 
formulating  a  truly  dynamic  version  of  the  Stratford 
criterion,  the  results  presented  in  this  manuscript  are 
mostly  limited  to  its  quasi-steady  implementation. 
Typical  separation  lines  obtained  in  this  mariner  are 
provided  in  Fig.  2.  For  the  body  geometry  described  in 
Section  6,  this  figure  shows  the  temporal  variations  of 
the  computed  separation  line  (i.e.,  the  separation  roll 
angle  .©s  measured  from  the  windward  stagnation  point 
as  a  function  of  streamwise  distance  along  the  body),  for 
a  step  change  in  angle  of  attack  at  t  =  0,  from  a  =  15°  to 
a  =  20*. 

In  each  crossflow  plane,  the  strength  of  the 
vorticity  jet  is  derived  from  the  observation  that  for  a  flat 
plate  boundary  layer  of  thickness  6,  the  streamwise 
vorticity  flux  per  unit  span  is  given  to  the  order  of 
boundary  layer  theory  by: 

*  o2 

]  undy  -  ~  (15) 
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where  Ue  is  the  velocity  at  the  edge  of  the  boundary 
layer.  At  the  point  of  separation,  it  is  assumed  that  a  net 
fraction  X  of  this  vorticity  flux  is  injected  into  the 
freestream.  This  method  has  been  formally  shown1  to  be 
equivalent  at  subsonic  speeds  to  the  "vortex  cloud" 
method  used  by  Mendenhall  and  Perkins.9  In  its 
numerical  implementation,  the  incremental  normal 
vorticity  flux  released  between  x  and  x  +  Ax  is  derived 
from  the  condition: 

x+Ax 

|  j  Vnnd3  '  §(V«p+WapUAx  (16> 

x  ax 

where  the  subscript  "sp”  denotes  a  value  at  separation, 
and  Sj  and  s2  are  values  of  the  curvilinear  coordinate 
along  the  body  in  the  crossflow  plane,  placed  on  either 
side  of  the  separation  point.  The  empirical  "vortex 
reduction  factor*'  X  determines  the  amount  of  vorticity 
shed  at  that  point  and  was  set  to  be  equal  to  0.6  in  the 
present  calculations. 

The  transfer  of  the  vorticity  jet  boundary  condition 
to  the  computational  fuselage  and  the  transport  of 
vorticity  from  that  point  are  analyzed  next.  For 


consistency  with  the  model  used  to  transfer  velocity 
boundary  conditions,  the  StrAtford  criterion  is 
implemented  on  the  true  surface  of  the  body  by  making 
use  of  the  inverse  transformation  Z  «  x' *(<*)  to  obtain  Cp 
(see  Section  5).  This  same  transformation  can,  in  turn,  fe 
used  to  compute  the  intersection  of  the  stagnation 
streamline  with  the  computational  boundary.  This 
location  would  determine  the  placement  of  the  vorticity 
jet  on  the  computational  fuselage.  At  tha  present  time, 
the  computer  code  is  set  up  to  prescribe  any  distribution 
of  normal  vorticity  fluxes  along  the  computational 
fuselage  boundary.  For  simplicity,  the  results  presented 
in  this  paper  consider  only  the  model  case  of  a 
concentrated  vorticity  flux  placed  at  the  closest  grid 
point  from  the  location  of  separation. 

An  additional  difficulty  in  the  integration  of  the 
vorticity  transport  equation  arises  from  the  fact  that,  in 
the  original  CAP-TSD  code,  vorticity  was  initially  pulled 
away  from  the  computational  domain  and  towards  the 
body,  due  to  the  downwash  effect.  This  problem  was 
circumvented  by  placing  the  body  parallel  to  the  grid  and 
considering  the  freestream  flow  to  be  inclined  at  an  angle 
of  attack  a.  As  a  first  approximation,  the  disturbance 
potential  0  can  be  solved  using  the  method  of  integration 
in  CAP-TSD,  as  long  as  a  remains  small,  although  in  the 
present  study  the  theory  is  pushed  well  beyond  its 
limits,  i.e.,  for  finite  a.  The  potential,  0  is  therefore 
replaced  by  (0  +  sin(a)  z)  in  order  to  resolve  the 
transport  of  vorticity  away  from  the  computational 
fuselage. 

b)  Farfield  and  Symmetry  Plane .  The  farfield 
boundary  conditions  are  based  on  characteristic  'non¬ 
reflecting'  boundary  conditions  for  0,  and  on  the  fact 
that  the  rotational  flow  component  must  vanish  away 
from  the  solid  boundary.3  The  non-reflecting  boundary 
conditions  are  derived  after  Whitlow.10  These  are: 

Upstream  boundary:  0-0 

Downstream  boundary :0X  -  -  jjsJ^t 

Upper  boundary:  0Z  -  0t  +  ^y  '(17) 

Lower  boundary:  0Z  -  2  +  Ay 

Spanwise  boundary:  0y  -  0t  -  Az 

where: 

A  -  M2  B  -  2M2  C  -  E+2F0 

eo  «>  X 

2  (18) 

D  -  <|ma>1/2  E  -  l-M2  P  -  -|(y+l)M2 

The  present  study  is  restricted  to  pitching  motion  only 
About  the  symmetry  plane,  v  and  n  are  antisymmetric  in 
y  therefore  at  y  =  0: 

0y  +  Az  =  0 ;  O  =  0  (19) 

c)  Horizontal  Fins  and  Wake .  For  horizontal 
lifting  surfaces,  the  boundary  conditions  on  0  are 
unchanged  from  their  original  treatment  in  CAP-TSD 
(Refs.  11  and  12).  Along  the  top  (+)  and  bottom  (-) 
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surfaces  of  the  surface,  the  vertical  disturbance  velocity 
is  specified  according  to  a  flow  tangency  condition: 

**  =  6±(x,y,t)  (20) 

where  6  represents  the  local  vertical  deflection  angle  at 
the  fin  surface.  Across  the  wake,  the  vertical  velocity 
component  $7  is  required  to  be  continuous,  while  the 
wake  circulation,  r,  is  convected  downstream12  by 
integrating  the  equation 

rx  +  rt=o  (21) 

It  may  be  shown13'14  from  a  contour  integral  in  the 
trailing  edge  plane  of  the  fin  that  an  initial  condition  at 
the  trailing  edge  is  given  by: 


rtc(y,0  =  rte(y,t)-^te(y,t)  (22) 

Equation  21  can  be  used  to  represent  at  each  streamwise 
location  the  effect  of  the  wake  on  the  body,  by  assigning 
a  strength  I*,  to  point  vortices  placed  at  each  spanwise 
grid  location  of  the  wake  cut.  This  horizontal  array  of 
point  vortices  is  then  mapped  into  o-space,  inducing  a 
complex  velocity 

§  "  1  25f  M*,-  <23> 


In  order  to  obtain  Cp  at  the  body  surface,  it  is 
necessary  to  know  q2  and  df  / dt  at  the  body.  This  is 
achieved  by  inverting  the  transformation  used  to  transfer 
boundary  conditions  from  the  true  body  surface  to  the 
cartesian  computational  fuselage.  Grid  points  Z 
belonging  to  the  computational  boundary  are  first 
projected  orthogonally  onto  points  Zt  which  belong  to 
the  true  body  surface.  Based  on  the  analytic  solution  for 
a  point  vortex  in  a  crossflow,  the  complex  velocities  at 
the  corresponding  points  in  <r-space  can  be  related 
simply.  It  can  be  shown  that  the  velocities  at  a  and  <r$ 
differ  by  an  additive  complex  quantity  K  which  depends 
only  upon  <r,  (T& ,  and  tne  strength  and  position  of 
potential  vortices  present  in  the  model.  Following  the 
derivations  of  Ref.  3,  this  same  transfer  function  K,  is 
used  to  obtain  the  complex  velocity  at  <r#,  based  on  the 
actual  (i.e.,  calculated)  velocity  at  a,  i.e.: 


dF 

dcr 


-1, 


<26> 


+  K(<T, CTs,<ro, or.), ro,rj> 

where  <r  is  the  mapped  location  of  the  rotational  body 
vortex  of  strength  ro,  and  the  <r/s  and  T/s  represent  the 
location  and  strength  of  the  w :  wakewortices.  The 
complex  transfer  function,  k,  is  given  by: 


where  <r :  is  the  location  of  the  individually  mapped 
vortices  and  ZT.  denotes  the  complex  conjugated  <7-.  The 
summation  is1  carried  out  over  the  set  of  spanwise  grid 
locations  representing  the  wake.  The  complex  velocity 
dF'/dff  induced  by  the  wake  circulation  is  subsequently 
added  to  dF/do-  in  Equations  12  and  13  to  give  the 
proper  spanwise  and  normal  boundary  conditions  on  the 
computational  fuselage. 


^aF^j  +  a~al 

J  (27) 

_ U] 

The  scaled  crossflow  velocity  at  the  body  surface  is 
given,  therefore,  by  the  inverse  mapping: 


5.  Pressure  Coefficient  Calculation 


At  the  solid  surfaces  the  pressure  coefficient  is 
calculated  using  the  isentropic15  relation: 


-w  -iv 
s  s 


(28) 
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where  q2  is  the  squared  modulus  of  the  total  velocity. 
Nixon  et  al.2  pointed  out  that  in  order  to  get  the  correct 
effect  resulting  from  the  tilting  of  the  lee-side  vortex 
away  from  the  surface,  it  is  necessary  to  apply  a  second 
order  correction  to  q2  so  that  it  includes  the  rotational 
component  of  streamwise  velocity: 

-  j  ih dz  (25> 


Similarly,  the  complex  potentials  at  <7g  and  a  can 
easily  be  related  using  the  following  relation: 


-  Sin«x>2|%e(<rs-<T))  J 


Hence,  a 0/at  I  s  is  not  equal  to  30/at.,  and  must  be 
modified  in  order  to  account  for  the  temporal  variations 
of  <70,.ro  and  T. ,  which  are  the  source  of  additional 
phase  lags'in  the  determination  of  the  pressure 
coefficient. 


The  inclusion  of  this  higher  order  term  has  been 
shown2'3  to  yield  more  accurate  results  and  lays  a 
theoretical  foundation  for  the  "empirical  correction"  used 
by  Mendenhall  and  Perkins.9 


6.  Results 

Time-accurate  calculations  were  performed  for  a 
missile-shaped  body  of  11:1  aspect  ratio  and  circular 
cross-section.  The  body  is  an  ogive-cylinder  with  a  3- 
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caliber  nose  and  8-caliber  straight  section,  chosen  to 
allow  comparisons  with  Tinling  and  Allen's16 
experimental  data.  In  order  to  demonstrate  the  current 
capabilities  of  the  code,  generic  rectangular  horizontal 
and  vertical  fins  could  also  be  added  to  the  body,  as 
shown  in  the  perspective  view  of  Fig.  3.  Each  fin  was 
infinitely  thin  ana  placed  at  zero  angle  of  attack  with 
respect  to  the  body.  The  leading  ana  trailing  edges  of 
the  fins  were  located  at  Xj*e/L  «  0.74  and  Xt  /L  =  0.86, 
respectively.  The  horizontal  lifting  surface -had  a  span 
equal  to  85%  of  the  maximum  body  radius.  The  body 
was  configured  at  a  mean  angle  of  attack  a  =  15*,  which 
could  be  oscillated  sinusoidally  with  an  amplitude  of 
+20%.  The  Mach  number  was  varied  between  Mw  =  0.3 
andM^sl.5. 

For  all  of  the  results  presented  in  this  manuscript, 
the  computational  mesh  size  was  45  x  32  x  50  in  the 
streamwise,  spanwise,  and  normal  directions.  The 
cartesian  grid  was  designed  to  optimize  spatial 
resolution  around  the  separated  flow  region,  and 
extended  at  least-five  body  diameters  in  the  spanwise 
and  normal  directions  (Fig.  1). 

In  order  to  verify  the  accuracy  of  the  analytic 
transfer  of  boundary  conditions  described  in  Section  4, 
pressure  distributions  were  computed  for  a  purely 
irrotational  calculation.  Fig.  4  shows  the  Cp  distribution 
as  a  function  of  roll  angle,  in  a  streamwise  plane  which 
corresponds  to  the  straight  section  of  the  body.  The  data 
points  (open  symbols)  represent  the  pressure  coefficients 
calculated  according  to  the  method  described  in 
Section  5,  where  the  grid  points  closest  to  the  actual 
body  (i.e.,  the  inner  ’‘corners”  in  Fig.  1)  have  been  used. 
This  computed  pressure  distribution  is  compared  to  the 
analytical  solution  for  potential  flow  around  a  circular 
cylinder  (dashed  line  in  Fig.  4).  Since  the  body  diameter 
is  constant  around  this  streamwise  location,  the  two 
pressure  distributions  are  expected  to  compare  relatively 
well,  provided  that  the  appropriate  vertical  shift  is 
added  to  the  analytical  solution.  Note  that  this  shift 
must  be  introduced,  in  order  to  take  into  account  the 
local  value  of  a  streamwise  component  of  velocity,  which 
differs  from  unity.  As  may  be  seen  from  Fig.  4,  the 
agreement,  in  the  absence  of  vorticity,  is  quite 
satisfactory.  This  implies  that  the  transfer  of  boundary 
conditions  that  was  derived  in  Ref.  3  may  be  used 
reliably,  to  the  extent  that  the  surface  rotational 
velocities  induced  by  the  distributed  vorticity  can  be 
adequately  represented  by  the  effect  of  a  point  vortex. 

The  decomposition  of  the  flowfield  into  rotational 
and  irrotational  velocity  components  is  illustrated  in  Fig. 
5,  for  a  crossflow  plane  intersecting  the  fins.  The  left 
hand  graph  shows  the  irrotational  velocity  component, 
+  sin(a)k,  where  k  is  the  unit  vector  in  the  positive  z- 
direction.  The  rotational  component  VxA  (middle  graph) 
illustrates  the  presence  of  a  formed  vortex  away  from  the 
surface.  The  total  velocity  field*(right~hand  side)  is  the 
superposition  of  rotational  and  irrotational  flow 
components. 

Fig.  6  shows  the  details  of  vorticity  contours  at 
X/D  =  10.1  in  the  absence  of  fins,  and  illustrates  the 
entrainment  of  vorticity  from  the  point  of  shear  layer 
separation,  towards  the  core  of  the  vortical  structure 
exhibited  in  Fig.  5.  The  corresponding  global  dynamics 
of  three-dimensional  vorticity  transport  and  vortex  roll¬ 
up  are  illustrated  in  Fig.  7  for  the  case  of  a  step  change  in 
angle  of  attack  at  t  =  0,  from  a  =  15®  to  a  =  20  •,  at  a 


Mach  number  M„  =  0.3.  In  the  present  calculations,  the 
incorporation  of  dynamic  effects  through  the  coupling  of 
outer  flow  dynamics  with  the  instantaneous  motion  of 
the  separation  line  was  inhibited  along  the  forebody.  In 
the  absence  of  this  restriction,  vorticity  at  the  nose  is 
introduced  too  far  away  from  the  surface.  This  results  in 
loss  of  coupling  and  early  separation  of  a  weak  tip 
vortex.  The  source  of  this  problem  is  similar  to  that 
identified  earlier  in  the  crossflow  plane  and  has  to  do 
with  the  cylindrical  nature  of  the  computational  fuselage 
in  CAP-TSD.  Consequently,  vorticity  jets  were  only 
introduced  from  the  end  of,  the  forcbody  section  (i.e., 
beyond  X/D  =  2.5),  as  seen  from  Fig.  7. 

The  anatomy  of  the  flowfield  is  perhaps  best 
described  by  the  evolution  of  the  velocity  field  at  various 
downstream  positions  along  the  body  (Fig.  8).  At 
X/D  =  4.0,  little  vorticity  has  been  injected  into  the  outer 
flow.  The  formation  of  a  small  region  of  reversed  flow 
first  becomes  evident  around  X/D  =  5.6.  The 
phenomena  of  vortex  roll-up  and  strengthening  with 
increasing  downstream  distance  can  be  observed  from 
the  evolution  of  the  flowfield  at  X/£)  =  7.9  and 
X/D  =  9.7.  These  observations  remain  qualitatively 
similar  to  those  of  Ref.  3,  although  in  the  latter  work,  the 
vortex  strength  and  its  vertical  location  were 
exaggerated,  due  to  the  inaccuracies  involved  in 
positioning  a  concentrated  vorticity  jet  at  a  distant 
computational  boundary;  as  previously  discussed. 


A  quantitative  comparison  of  pressure 
distributions  at  several  downstream  locations  is  made  in 
Fig.  9.  These  pressure  distributions  are  compared  to  the 
experimental  data  of  Tinling  and  Allen16  (open  symbols) 
and  the  predictions  of  Mendenhall  and  Lesieutre1 7 
(dashed  lines),  using  the  "vortex  cloud’*  method.8  At  all 
three  locations,  the  current  predictions  are  seen  to 
overshoot  both  the  experimentally  acquired  values  and 
the  vortex  cloud  predictions  in  the  separated  flow 
region.  By  comparison  with  the  results  of  Nixon  et  al.,3 
one  can  infer  that  the  discrepancy  may  be  due  to  the 
lower  strength  of  the  rotational  flow  component.  In  the 
vortex  cloud  method,  incipient  separation  was  reported 
to  take  place  around  X/D  » 1.0,  out  was  only  allowed 
beyond  X/D  a  2.5  in  the  present  calculations.  The  level 
of  disagreement  between  the  experimental  data  and  the 
current  results  is  also  seen  to  diminish  with  increasing 
downstream  distance  (a  result  also  observed  in  Ref.  2). 


Fig.  10  illustrates  the  effect  of  Mach  number  at  an 
angle  of  attack  a  -  15*  by  comparing  crossflow 
velocities  /  vorticity  distributions  at  Mw  =  0.3  and 
M„  =  15.  As  can  be  seen  in  this  figure,  the  vorticity  shed 
into  the  lee-side  vortex  is  weaker  in  the  supersonic  case. 
This  results  in  a  smaller  vortex,  residing  closer  to  the 
surface  than  in  the  case  of  subsonic  flow. 


One  of  the  objectives  of  this  study  is  to  assess  the 
aerodynamic  phase  lags-which  occur  in  unsteady 
separated  flow.  The  phase  lags  result  from  various 
dynamic  effects,  which  include  Tags  in  the  position  and 
strength  ofyortices  with  respect  to  the  body  motion. 
These  phase  lags  originate  from  the  fact  that  the  position 
and  strength  of  a  vortex  at  any  streamwise  location 
along  the  body  are  the  integral  consequence  of  vorticity 
transport  dynamics  on  the  one  hand,  and  the  unsteady 
rate  at  which  vorticity  is  fed  into  the  vortex  sheet  on  the 
other.  Therefore  the  'history  effect'  is  related  to 
convective  phase  delays,  as  well  as  additional  time  lags 
associated  with  the  separation  process  itself.  Ref.  8 


-6- 


provides  a  detailed  account  of  the  application  of  indicial 
theory  [ o  the  prediction  of  unsteady  separation. 

In  this  approach,  the  motion  of  the  separation 
point  and  the  unsteady  vorticity  flux  in  a  given  crossflow 
plane  are  analyzed  using  two-dimensional  time-accurate 
Navier-Stokes  calculations.  The  study  described  in  Ref.  8 
shows  that  for  viscous  flow  about  a  two-dimensional 
cylinder,  the  time-dependent  behavior  of  both  separation 
angle  and  vorticity  flux  can  be  accurately  predicted 
using  a  convolution  integral,  based  on  the  knowledge  of 
the  step  response  of  the  flowfield  to  a  smalt  perturbation 
in  Macn  or  Reynolds  number. 


One  of  the  results  was  that  periodic  fluctuations  of 
the  vorticity  flux  exhibited  a  50%  amplification  at  a 
reduced  frequency  based  on  diameter  and  crossflow 
velocity  of:  kf  =  2rtfD/Uc  =  23.  For  the  11:1  aspect  ratio 
body  described  above  and  at  a  mean  angle  of  attack 
a  =  15",  this  corresponds  to  a  reduced  frequency  based 
on  body  length  of  approximately  k  =  6.5.  At  a  Mach 
number  M,,,  =  0.9  ana  for  a  hypothetical  body  length 
L  =  10  m,  this  would  correspond  to  a  frequency 
/  =  10.5  Hz,  which  is  well  within  the  range  of  structural 
frequencies  exhibited  by,  e.g.,  soa-skimming  missiles.18 


To  illustrate  this  effect,  the  total  normal  force 
coefficient  (i.e.,  including  body  and  fins)  was  recorded 
(Fig.  11)  as  a  function  of  time,  for  the  case  of  an 
oscillating  angle  of  attack  a(t)  =  a0  +  ocjHfOsinCktU^/L), 
where  H(t)  is  the  Heaviside  step  function,  aQ  =15*, 
ctj  =  3*,  and  k  =  6.47.  Two  calculations  were  performed 
at  a  Mach  number  =  0.9.  In  the  first  calculation  (solid 
line),  the  Stratford  criterion  (Equation  14)  is  used  at  each 
time  step  to  determine  the  separation  location.  This 
location  is  used  to  obtain  V,p  and  Wsp  in  Eq.  16.  The 
value  of  the  crossflow  velocity  at  Separation  then 
specifies  the  normal  vorticity  flux  according  to 
Equation  16.  This  first  scenario  represents,  t  before,  a 
"quasi-stead/  implementation  of  the  Stratford  criterion. 


In  a  second  calculation  initiated  from  the  same 
steady  state,  a  negligible  amplification  of  separation 
angle  fluctuations  was  assumed,  while  a  time-dependent 
vorticity  flux  y(x,t)  was  imposed  based  onthe  findings  of 
Ref.  8.  This  flux  had  an  identical  mean,  y(x),  to  that  of 
the  first  calculation,  but  had  a  fluctuating  component 
equal  to  /5(y(x, t)*y(x)),  where  p  =  135.  The  amplification 
factor  p  was  chosen  to  be  real  because  previous  results8 
have  indicated  that  the  phase  delay  at  that  frequency  is 
close  to  zero.  The  resulting  normal  force  coefficient  is 
indicated  by  the  dashed  line  in  Fig.  11.  As  may  be  seen 
from  this  figure,  the  total  normal  force  lags  the  motion 
by  approximately  320«  in  both  cases.  However, 
significant  differences  are  observed  between  peak  values 
and  harmonic  content  of  the  oscillatory  loads. 


As  a  final  word  of  caution,  it  should  be  noted  that 
the  results  of  Fig.  11  do  not  represent  a  realistic 
evaluation  of  the  importance  of  dynamic  effects  on  the 
separation  process,  but  are  simply  meant  as  an  example 
illustrating  the  way  in  which  small  changes  in  the 
separation  (vorticity)  boundary  conditions  may  affect  the 
total  loads.  Among  the  reasons  that  this  second 
calculation  may  not  be  realistic  are  the  fact  that  the 
results  of  Ref.  8  are  limited  to  low  Reynolds  numbers 
and  small  amplitudes.  Additionally,  the  instantaneous 
perturbations  of  the  Navier-Stokes  equations  involved 
step  changes  in  Mach  number  and  Reynolds  number, 
rather  than  in  the  surface  pressure  distribution,  which  is 


the  basis  of  the  Stratford  criterion.  In  other  words,  the 
quasi-steady  implementation  of  the  Stratford  criterion 
does  not  correspond  to  a  quasi  steady  change  in  angle  of 
attack,  as  evidenced  by  the  temporal  evolution  seen  in 
Fig.  2.  In  spite  of  these  differences,  it  should  be  noted 
that  additional  motion  of  the  separation  point  has  been 
ignored,  and  that  a  seventeen-fold  increase  in  Reynolds 
number  (Fig.  12)  suggests  an  increase  in  phase  delays,  as 
well  as  further  amplification  of  fluctuations  in  the 
vorticity  flux.  This  implies  that  the  differences  in  the 
predictions  of  Fig.  11  may  actually  underestimate  the 
effects  of  unsteady  separation. 

7.  Concluding  Remarks 


The  theory  that  was  developed  in  Ref.  3  to  treat 
flow  separation  and  related  vortex  effects  in  unsteady 
transonic  flow  around  slender  bodies  was  implemented 
using  a  modified  version  of  the  CAP-TSD5  computer 
code.  This  theory  involves  the  simultaneous  solution  of 
a  modified  Transonic  Small  Disturbance  equation,  a 
vector  potential  equation,  and  a  three-dimensional 
unsteady  vorticity  transport  equation3.  In  the  present 
work,  refinements  in  the  representation  of  the 
computational  fuselage  are  shown  to  yield  significant 
improvements  over  previous  predictions.  The  results  of 
time-accurate  computations  for  complete  missile 
configurations  at  subsonic,  transonic,  and  supersonic 
speeds  suggest  that  realistic  angle  of  attack 
configurations  may  be  calculated  using  CAP-TSD,  thus 
showing  considerable  potential  for  aeroelastic 
computations  and  unsteady  aerodynamics. 
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Figure  2.  Time-Dependent  Separation  Line  for  a  Step 
Change  in  Angle  of  Attack  at  t  =  0,  from  a  =  15-  to 
a  =  20»(M  =0.3). 


Figure  3.  Perspective  View  of  Prototypical  Slender  Body 
Under  Consideration,  Including  Horizontal  and 
Vertical  Fins. 
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Figure  1.  Cross-sectional  Views  of  the  Cartesian  Grid 
(Left:  Full  Grid ;  Right:  Close-up  Illustrating 
Serrated-Edged  Computational  Fuselage, 
with  Roll  Angle  Convention). 
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Figure  4.  Illustration  of  the  Accuracy  of  the  Transfer  of 
Boundary  Conditions  by  Comparison  of  Pressure 
Coefficient  Distributions  between  a  Purely 
Irrotational  Calculation  (X/D  =  7.9)  and  the 
Analytical  Two-Dimensional  Potential  Flow 
Solution  ( a  =  15*,  M^  =  0.3). 
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Figure  5.  Illustration  of  the  Helmholz  Decomposition  Procedure  at  the  Leading  Edge  of  a 
Horizontal  Lifting  Surface  (X/L  =  7.4,  =  0.3,  a  =  15*). 
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Figure  6.  Vorticity  Distribution  at  X/D  =  10.1  (a  =  15°, 
M  =0.3). 


Figure  7.  Half-Space  Representation  of  Streamwise 
Vorticity  Contours  Illustrating  the  Time- 
Dependent  Formation  of  the  Leeward  Side 
Vortical  Structure,  for  a  Step  Change  in  Angle  of 
Attack  at  t  =  0,  from  a  -  15*  to  a  ~  20* 
(M=0.3). 
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Figure  8. 
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Streamwise  Evolution  of  the  Vortical  Crossflow  (M^ 
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Figure  12.  Reynolds  Number  Comparison  of  Theoretical  Predictions  (Based  on  Indicial  Theory)  of 
Phase  Delays  (top  graph)  and  Amplification  Factors  (bottom  graph)  of  the  Vorticity  Flux, 
Subject  to  Cross-Flow  Mach  Number  Oscillation  About  Mc  =  0.25.  Open  Symbols  Are  the 
Result  of  Numerical  Experiments  Based  on  Two-Dimensional  Navier-Stokes  Calculations.8 
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